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ABSTlUkCT 


Application  of  classic  triangulation  methods  will  allow 
the  location  of  a  radar  to  be  determined  by  passive  sensors. 
Through  the  use  of  modern  digital  signal  processing  techniques 
this  estimate  can  be  made  in  a  simpler  fashion  using  a 
conventional  receiver. 

In  this  thesis  a  technique  is  developed  for  time 
difference  of  arrival  (TDOA)  estimation  using  a  frequency 
domain  based  correlation  detector  driven  by  an  envelope 
detector.  Time  lag  boundaries  are  defined  on  the  output  of  the 
correlator.  A  fixed  detection  threshold  is  calculated  to 
permit  constant  false  alarm  rate  (CFAR)  detection.  The 
performance  of  the  correlation  detector  is  plotted  as  a 
receiver  operating  characteristic  (ROC)  curve  as  a  function  of 
signal  to  noise  ratio  (SNR) .  An  interactive  MATLAB  software 
program  is  provided  to  perform  either  spectral  domain  or  time 
domain  based  correlation. 

Spectral  domain  based  correlation  uses  the  Fast  Fourier 
Transform  (FFT) .  Implicit  with  the  use  of  the  FFT  are  finite 
arithmetic  internal  processing  errors  which  are  modeled  as 
independent  uncorrelated  noise  sources.  A  method  is  presented 
to  account  for  SNR  degradation  at  the  output  of  the  FFT. 
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I .  INTRODUCTION 


A.  BACKGROUND 

Interception  of  radar  transmissions  by  passive  sensors  is 
one  of  the  responsibilities  of  electronic  intelligence. 
Analysis  of  the  signal  of  a  radar  can  provide  an  estimate  of 
the  distance  and  direction  of  an  emitter  relative  to  an 
intercepting  receiver.  Knowing  the  location  is  an  important 
consideration  in  the  evaluation  of  strategic  and  tactical 
cap2Q3ilities  of  the  radar. 

Triangulation  is  a  traditional  technique  where  the  angle 
of  arrival  of  a  radar  signal  at  two  spatially  separated 
receivers  provides  an  estimate  of  the  transmitters  location. 
With  the  advent  of  the  computer,  digital  signal  processing 
algorithms  can  be  implemented  that  provide  emitter  range 
information  and  hence  allow  localization  in  a  sequential 
sense.  One  of  these  processing  algorithms  is  called  the  time 
difference  of  arrival  (TDOA)  method. 

B.  TIME  DirnERSNd  OW  ARRIVAL 

The  TDOA  algorithm  relies  on  two  bistatic  passive  sensors 
attempting  to  receive  the  transmission  of  an  active  radar. 
Reception  by  the  two  sensors  is  preferentially  done  in  the 
main  lobe  of  the  transmitting  radar,  but  due  to  the 
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nonsynchronous  nature  of  the  task,  reception  through  one  of 
the  sidelobes  of  the  transmitting  radar  is  expected  with  a 
resulting  decrease  in  received  signal  power. 

Given  a  transmitted  radar  pulse  f (t) ,  two  physically 
separated  passive  sensors  will  roost  likely  receive  this  pulse 
at  times  ti  and  ti+T.  The  time  difference  T  is  proportional  to 
their  differential  distance  relative  to  the  emitter.  This  time 
difference  of  arrival  can  be  used  in  the  estimation  of  the 
differential  range  of  the  emitter  to  the  receivers.  The 
reception  of  the  first  pulse  in  a  pulse  burst  by  the  sensors 
is  paramount  in  performing  the  TDOA  estimate. 

The  concept  of  measuring  the  TDOA  of  a  radar  pulse  between 
two  separated  receivers  can  be  modeled  by  a  hyperbola.  The  two 
sensors  whose  positions  are  known  are  located  at  the  two  focal 
points  of  the  hyperbola.  A  transmitting  radar  is  located  in 
the  area  surrotinding  the  two  sensors.  Figure  1  shows  this 
model . 

A  hyperbola  is  the  locus  of  points  in  which  the  distance 
from  euiy  one  point  on  the  hyperbola  to  one  focus  differs  by  a 
constant  amount  from  the  distance  to  the  other  focus.  This 
differential  distance  is  proportional  to  a  difference  in 
arrival  time  t,  ol  a  radar  pulse  as  detected  between  the  two 
sensors.  Once  t  is  measured  the  hyperbola  can  be  drawn.  If  the 
sensors  (moving  aircraft  or  nongeostationary  satellites) 
repeat  the  differential  arrival  time  measurements  at  several 


2 


time  intervale,  a  series  of  hyperbolas  will  be  generated, 
fixing  the  location  of  the  transmitting  radar.  Alternately, 
the  angle  of  arrival  at  either  sensor  can  be  used  to  define 
the  position  of  the  transmitter  on  the  hyperbola.  A  good 
introduction  into  the  use  of  the  hyperbola  in  navigation  is 
found  in  [Ref.  l:p.  27]. 


Range  (mete'Si 


rigure  1.  Time  difference  of  arrival  model. 

This  discussion  assumes  that  the  radar  signal  received  by 
the  two  sensors  is  induced  by  the  direct  wave  from  the 
transmitter.  The  effects  of  refraction  and  reflection  on 
propagating  radar  waves  are  not  covered  here. 
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The  purpose  of  this  thesis  is  to  develop  and  test  an 
algorithm  to  estimate  the  differential  arrival  time  T,  of  a 
pulsed  radar  signal  collected  by  two  passive  sensors.  If  the 
pulse  burst  is  collected  in  the  time  domain,  time  correlation 
can  be  used  to  generate  a  TDOA  estimate.  The  Fast  Fourier 
Transform  (FFT)  can  also  be  used  to  perform  the  correlation  of 
the  two  received  signals  (i.e.,  fast  correlation). 

Additive  noise  n(t)  ,  superimposed  onto  a  transmitted  radar 
signal  f (t)  by  the  environment  and  the  radar  receiver,  creates 
a  composite  signal  s (t)  at  the  receiver 

s  ( t)  =  f { t)  +  n  ( t)  .  (1) 

The  use  of  the  Fourier  Transform  for  TDOA  estimation  is 
attractive  because  of  the  processing  gain  (P6)  a  radar  signal 
receives  in  being  transformed  from  the  time  domain  to  the 
frequency  domain  during  the  detection  process.  FFT  processing 
gain  for  real  valued  signals  can  be  approximated  by 

PG  (dB)  =  (logj  [number  points  in  data  record]  -i)  *3  .  (2) 

See  for  example  [Ref.  2:p.  33].  The  signal  to  noise  ratio 
(SNR)  at  the  output  of  a  FFT  (SNR<^)  is  a  function  of  the  PG 
and  the  SNR  at  the  input  to  the  FFT  (SNRj„)  .  This  relation  is 
given  by 
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(3) 


SNRo^  =  SNR,^  -  PG  . 

The  corruptive  effects  of  additive  noise  are  reduced  by  the 
transformation,  allowing  the  frequency  domain  detection  of  the 
signal . 

C .  RADAR  TYPES 

The  TDOA  algorithm  is  primarily  applied  to  pulsed  radars. 
The  continuous  wave  (CW)  radar  provides  target  radial  velocity 
through  the  Doppler  shift  of  its  return.  The  CW  radar  does  not 
use  a  pulsed  transmission.  Therefore,  the  TDOA  algorithm 
cannot  not  be  directly  applied  to  the  reception  of  CW  signals . 

Pulsed  radars  can  be  divided  into  two  broad  categories, 
ordinary  low  pulse  repetition  frequency  (PRF)  pulsed  radars 
and  the  pulse  Doppler  radars. 

Low  PRF  radars  can  yield  unambiguous  range  information, 
while  their  radar  returns  do  not  give  target  velocity 
information  directly.  Generally  because  of  the  long  duty  cycle 
of  the  low  PRF  radar,  the  probability  of  receiving  all  radar 
pulses  in  a  pulse  burst  by  two  passive  sensors  is  great.  The 
TDOA  algorithm  can  then  estimate  the  differential  distance  to 
the  emitter  using  the  low  PRF  burst. 

The  pulse  Doppler  radar  uses  coherent  transmission  and 
reception  with  a  moderately  high  PRF  [Ref.  3:p.  17.1].  it 
gives  both  rauige  and  radial  velocity  information  about  a 
target,  though  not  with  the  seune  accuracy  of  the  ordinary 
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pulsed  radar  or  CW  radar  respectively.  In  the  pulse  Doppler 
radar,  as  in  all  non-bistatic  radars,  the  receiver  must  be 
turned  off  while  the  transmitter  is  on.  Because  of  this,  a 
range  ambiguity  exists.  Echoes  having  delay  times  equal  to  an 
integral  multiple  of  the  PRF  will  be  undetected.  Target  echoes 
that  remain  within  this  blanked  time  interval  will  remain 
undetected.  A  second  source  of  range  eunbiguity  exists.  For  a 
target  located  at  a  distance  r^  from  the  tracking  radar,  all 
other  targets  in  the  main  beam  of  the  radar  a  distance 


2r^,  3r^,  4r^,  ....  kr^ 
where  k  =  integer  , 


(4) 


will  appear  to  be  approximately  the  aanie  distamce  from  the 
radar . 

Pulse  Doppler  radars  can  be  divided  into  a  low,  medixim  and 
high  PRF  class .  Both  the  low  emd  medium  PRF  pulse  Doppler 
radar  treuismissions  can  be  passively  intercepted  and  a 
relatively  simple  TDOA  estimate  cam  be  formed. 

Reception  of  each  pulse  in  a  pulse  burst  from  a  high  PRF 
pulse  Doppler  radar  by  two  independent  passive  sensors  cam  be 
difficult.  If  either  sensor  misses  the  first  pulse,  the  TDOA 
estimate  will  be  in  error.  This  is  true  for  any  radar,  but  has 
a  higher  likelihood  of  occurring  in  the  high  PRF  radars  due  to 
the  shorter  pulse  widths  used.  Should  the  pulse  burst  be  coded 
so  that  a  nonuniform  time  pulse  sequence  is  generated  (i.e.. 
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staggered  PRF) ,  then  missed  pulses  by  either  sensor  will  not 
seriously  degrade  the  validity  of  the  TDOA  estimate. 

D.  FAST  FOORISR  TRANSFORM  OUTPUT  SIOIAL  TO  HOISB  RATIO 

The  Fourier  transform  is  a  computationally  intensive 
mathematical  operation.  When  implemented  on  an  FFT  processor 
using  fixed  point  arithmetic,  internal  processing  errors  are 
generated.  These  errors  are  a  function  of  the  transform  size, 
and  the  number  of  bits  used  internally  in  the  FFT  processor. 

The  primary  building  block  of  the  FFT  is  the  butterfly 
algorithm.  There  are  three  dominant  processing  errors  that 
occur  within  the  butterfly.  They  are  scaling,  truncation  and 
trigonometric  errors.  These  can  be  modeled  as  independent, 
uncorrelated  noise  sources.  The  cumulative  effect  of  these 
noise  sources  is  to  reduce  the  potential  SNR  at  the  output  of 
the  butterfly,  and  therefore,  the  FFT  processor. 

This  thesis  will  also  examine  the  effects  of  noise  caused 
by  FFT  processing,  and  will  present  a  method  to  account  for 
the  degradation  of  the  output  SNR  due  to  finite  arithmetic. 
Certain  FFT  implementations  will  keep  the  signal  energy  larger 
with  respect  to  the  three  noise  power  sources  than  others . 
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II.  TIME  DOMA.IN  BASED  CORKELATION 


A.  BACKGROUND 

The  correlation  function  R(T)  measures  the  degree  of 
similarity  between  two  signals  s{t)  and  g(t)  and  is  described 
by 

R{t)  =  (t)g(-c  +  t)dt 
where  g(x  +  t)  =  g(c)  displaced  by  the  shift  t 

Correlation  can  be  performed  in  a  radar  receiver  between 
a  target  return  corrupted  by  additive  noise,  and  a  replica  of 
the  transmitted  signal  maintained  by  the  radar.  This  replica 
is  designed  into  the  frequency  response  of  the  matched  filter 
of  the  radar  receiver.  [Ref.  4:p.  373]. 

Correlating  two  signals  produces  an  output  that  does  not 
resemble  either  of  the  two  input  signals.  The  shape  of  the 
received  waveform  is  destroyed  during  the  correlation  process. 
The  useful  item  of  correlation  is  an  amplitude  peak,  where  the 
location  along  the  time  axis  of  this  peak  indicates  the  eunount 
of  time  that  one  signal  lags  the  other. 

If  the  correlation  process  occurs  between  a  signal  and  a 
delayed  replica  of  itself,  the  process  is  called 
autocorrelation.  If  two  dissimilar  signals  are  correlated,  the 
process  is  called  crosscorrelation. 
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B.  AUTOCORRBZATION 


Autocorrelation  is  the  process  by  which  two  identical 
signals,  one  a  delayed  replica  of  the  other,  are  correlated. 
The  autocorrelation  function  of  a  discrete  time  sequence  x(i) 
is  defined  by  [Ref.  5;p.  556] 

w-i-UI 

""  51  x(i)x(i  +  1) 

i-a  (6) 

where  1  =  shift  operator 

N  =  number  of  data  points  in  x in) 

This  function  is  not  scaled  with  respect  to  the  shift  operator 

nor  the  number  of  data  points.  Figure  2(a)  shows  the  effect  of 

autocorrelation  on  a  burst  of  13  pulses  of  unit  aunplitude. 

Normalized  autocorrelation  produces  a  triangular  envelope 

waveform  of  height  one  at  a  time  lag  of  zero.  The  location  of 

this  peak  indicates  that  the  time  shift  for  maximum  overlap 

between  the  pulse  burst  and  a  shifted  replica  of  itself  is 

zero. 

For  illustration,  a  Barker  coded  sequence  of  pulses  of 
length  13  is  also  autocorrelated  in  Figure  2 (b) .  The  Barker 
code  is  given  by  +++++ — ++-+-+,  where  +  and  -  could  represent 
0  and  It  radians  carrier  phase  shift,  respectively.  Barker 
coding  of  a  pulse  sequence  constructively  modifies  the 
correlation  output  to  magnify  the  time  lag  peak.  The  zero  lag 
point  has  the  greatest  aunplitude.  Adjacent  peaks  are 
significantly  reduced  in  amplitude.  Because  of  this  reduction. 


9 


Figur*  2.  (a)  Normalized  autocorrelation  of  13  equal 
aunplitude  pulses,  (b)  Normalized  autocorrelation  of  13  Barker 
coded  pulses. 

if  noise  is  added  to  the  received  signal,  a  Barker  coded  radar 
pulse  burst  will  have  a  correlation  peak  that  is  not  as  easily 
confused  with  adjacent  peaks. 

Now  consider  the  autocorrelation  of  a  random  process.  An 
independent,  identically  distributed  zero  mean  Gaussian 
sequence  with  variance  of  one  is  shown  in  Figure  3 .  The 
autocorrelation  of  this  sequence  produces  a  peak  at  zero  time 
lag.  There  is  no  strong  correlation  of  the  Gaussian  noise 
except  at  zero  time  lag.  The  autocorrelation  of  Rayleigh 
distributed  noise  is  shown  in  Figure  4 . 
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Rayleigh  distributed  noise  has  a  mean  (dc)  value  of 

mean  value  =  o 

where  a  =  standard  deviation 

The  autocorrelation  function  of  Rayleigh  noise  has  two 
components .  A  triangular  envelope  caused  by  the  dc  component 
of  the  Rayleigh  noise^  and  a  weighted  delta  function  at  the 
apex  of  the  triangle  which  indicates  the  time  shift  of  maximvim 
overlap  of  the  two  sequences.  Clearly  maximum  overlap  occurs 
at  a  time  shift  of  zero  for  any  autocorrelation  function. 

C .  CROSSCOBBBLATION 

Crosscorrelation  is  the  process  wherein  two  different 
signals  are  correlated.  The  crosscorrelation  function  of  a 
discrete  time  sequence  is  defined  by 

/f-i-UI 

^  x(i)y(i  +  1) 

where  xii)  and  y(i)  are  two  di  fferen t  sequences . 

This  function  is  not  scaled  with  respect  to  the  shift  operator 
nor  the  number  of  data  points.  Figure  5  shows  the 
crosscorrelation  of  two  uncorrelated  Gaussian  noise  sequences 
euid  two  uncorrelated  Rayleigh  noise  sequences.  Clearly,  no 
dominamt  amplitude  peak  is  visible  in  either  crosscorrelation 
plot.  Lack  of  a  dominant  peak  indicates  there  is  no  strong 
correlation  between  the  two  sets  of  noise  sequences. 
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Figure  5.  (a)  Normalized  Gaussian  noise  crosscorrelation, 

(b)  Normalized  Rayleigh  noise  crosscorrelation. 


D .  CORRKZATKMI  COBFriCIlNT 

A  measure  of  the  linear  dependence  between  any  two  zero 
mean  stationary  time  functions  x(t)  and  y(t)  is  the 
correlation  coefficient  function  and  is  defined  [Ref.  6;p.  48] 
as 

(  )  r^(T) 

(»> 

wheze  |Pxy(t)  |  i  1  for  all  x. 

This  function  is  the  crosscorrelation  function  normalized  by 
the  square  root  of  the  maximum  values  of  the  individual 
autocorrelation  functions.  The  normalizing  factor  is  not  a  lag 
dependent  quantity. 
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Given  two  time  series  x(t)  and  y(t)  where  a  is  any 
positive  number 

Pj^(t)  =  1  if  x{t)  =  a-y(c) 

if  =  -a*y(t)  (10) 

P^(t)  =  0  if  x{t)  and  y(t)  are  uncorrelated. 


Another  way  of  defining  a  normalized  correlation 
coefficient  for  a  discrete  time  series  of  finite  duration  is 
given  by 


P*v<-Z)  = 


^(-i)  y(i  +  1) 


i-0 


w-i 


(11) 


\  \  E 

N  N  i'O 


B.  SQOARB-LMV  OBTBCTIOM  AMD  CORBBLATIOM 

If  x(t)  is  a  real  valued  function  of  time,  the  Hilbert 
transform  of  x(t)  is  defined  by  [Ref.  6:p.  484]  as 


Hit)  =  H  [x(t)] 


x(u) 

X  ( t  ~  u) 


du 


=  x(t)  *  (^) 

X  t 

where  *  =  convolution  operation 


(12) 


The  output  of  a  square-law  envelope  detector  u(t),  can  be 
described  by 
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(13) 


u  ( t)  =  (c)  +  jc^  (  t) 

where  x(t)  =  real  valued  time  function 

input  to  square-law  detector 

The  correlation  coefficient  function  of  two  time  functions 
x(t)  and  y(t)  is  defined  in  Equation  9.  The  Hilbert  transform 
of  the  correlation  function  is  defined  as 

P^(x)  =  H  [p^(T)] 

O^Oy  (14) 

where  and  Oy  =  standard  deviation  of  x  and  y 

respectively  . 

The  correlation  coefficient  for  two  square-law  envelope 
detected  signals^  r  is  defined  by  [Ref.  6:p.  512]  as 


Puv(^)  *  ^ 

where  u{t)  *  square- law  detected  signal  of  x(t) 
v(  t)  *  square- law  detected  signal  of  y{t) 


(15) 


The  quantity  p„{X)  produces  a  sharper  correlation  peak  them 
p^(t)  .  This  sharpening  of  the  correlation  function  peak  output 
can  aid  in  locating  the  TDOA  correlation  peak.  An  example  of 
both  p^(t)  and  PaT('*)  is  plotted  in  Figure  6. 


r.  BHVBLOPK  COBBBZATIOH  AMD  BXZJITSD  STATISTICS 

In  this  thesis r  time  domain  based  correlation  will  be 
performed  on  the  output  of  an  envelope  detector.  An  envelope 
detector  can  be  easily  implemented  in  hardware  using  a  diode. 
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1 


normalized  correlation  coefficient 


1 


normalized  correlation  coefficient  and  Hilbert 


line)  . 

To  derive  the  expected  value  and  variamce  of  the 
crosscorrelation  function,  two  cases  must  be  considered.  These 
cases  are  noise  only  present,  and  signal  plus  noise  present  in 
both  channels  of  the  correlator.  Under  the  noise  only  case, 
the  statistics  for  the  correlator  output  are  easily  derived. 

We  consider  two  real,  independent  identically  distributed, 
zero  mean  noise  series  at  the  input  to  the  correlator.  Each 
series  has  zero  mean  and  variaiice  (T^  .  It  is  shown  in  ^pendix 
A,  that  the  expected  value  of  the  output  of  the 
crosscorrelation  fxanction  r,,  (f)  is 

Etr^(I)]  »  0  .  (16) 
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The  variance  of  the  output  of  the  crosscorrelation  function 
r,y(C)  is 


o2  =  (N-|i|)o^ 
where  N  =  number  of  data  points 


(17) 


Since  equal  variances  are  assumed  for  both  time  series,  c*  is 
the  square  of  the  variance  of  either  noise  series. 


G.  BLOCK  COBBBLKTIOM 

Two  time  sequences  can  be  correlated  by  any  of  three 
different  methods.  These  will  be  discussed  below. 

The  first  method  blocks  the  length  of  each  input  sequence. 
Correlation  of  these  blocks  will  provide  an  output  with  a 
variance  given  by  Equation  17.  The  linear  dependence  of  the 
correlation  variance  on  (N- | I | )  is  shown  in  Figure  5  for  zero 
mean  Gaussian  noise.  Figure  5  shows  the  correlated  output  at 

(zero  time  lag)  is  scaled  by  (N>|0|)»N.  Correlation  output 
at  a  time  lag  of  ±N  are  scaled  by  (N'>|n|)*0.  The  linear 
dependence  of  the  correlation  output  variance  on  the  time  lag 
does  not  allow  a  constant  detection  threshold,  which  is  needed 
to  automatically  determine  the  location  of  the  correlation 
peak. 

A  second  method  of  correlation  blocks  a  fixed  length 
sequence  x(n)  and  correlates  it  with  an  arbitrarily  long 
sequence  y (n) .  The  correlation  output  for  this  algorithm  has 
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a  variance  of  No*.  Clearly,  no  scaling  of  the  correlation 
output  variance  as  a  function  of  time  lag  would  occur. 

A  third  method  is  scaling  applied  to  the  output  of  the 
first  technique.  The  correlator  output  is  multiplied  by  a 
weighting  function,  l/(N-|Cj)^^*  to  correct  for  the  lag 
dependency  of  the  variance.  The  selection  of  a  constant 
detection  threshold  will  be  addressed  in  Chapter  V. 
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III.  FRSQUSNCY  DOMA.IM  BASED  CORRELATION 


A.  FOURIER  TRANSFORM 

Frequency  domain  techniques  can  be  used  to  perform  fast 
correlation  rather  than  the  time  domain  techniques  previously 
discussed.  The  Fourier  transform  translates  the  discrete  time 
sequence  x (n)  into  the  frequency  domain  . 

For  a  finite  duration  time  series  x(n),  the  discrete 
Fourier  transform  (DFT)  is  defined  by  the  pair  of  equations 
[Ref  5:p.  100] 

w-i  -ii2i)nk 

X(k)  -  ^x{n)  e  ''  .  0  i  k  N-1 

n’O 

(18) 

x{n)  =  ^  Xik)  e  ^  ,0  &  n  i  N-l  . 

^  jfso 

Equation  18  is  used  to  transform  the  N  discrete  seunple  values 
into  N  frequency  domain  coefficients. 

The  following  discussion  is  based  on  [Ref.  5:p.  286].  The 
DFT  is  a  computationally  intensive  algorithm  if  implemented 
directly.  The  direct  evaluation  of  the  DFT  requires  N*  complex 
multiplications.  The  time  required  for  the  evaluation  of  a  DFT 
is  proportional  to  N*. 

Algorithms  have  been  designed  that  take  advantage  of  the 
symmetry  and  periodicity  of  the  DFT,  reducing  the  eunount  of 
computation  required  to  solve  the  DFT.  Any  of  these  algoritlims 
are  referred  to  as  the  FFT. 
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Generally,  the  FFT  requires  (N/2)log2N  complex 
multiplications.  This  smaller  computational  workload  becomes 
significant  as  N,  the  niomber  of  data  points  in  the  time 
series,  becomes  larger. 

B.  NORMALIZED  FREQUENCY  DOMAIN  CORRELATION 

In  some  radar  receivers  the  processing  of  pulse  sequences 
occurs  in  the  frequency  domain.  Therefore,  a  frequency 
representation  of  the  time  domain  pulses  may  already  exist. 
The  frequency  coefficients  do  not  necessarily  need  to  be 
transformed  back  into  the  time  domain  to  perform  the 
crosscorrelation . 

Frequency  domain  data  can  be  crosscorrelated  by 
multiplying  one  set  of  coefficients  by  the  complex  conjugate 
of  the  second  set  of  coefficients.  The  product  is  translated 
back  to  the  time  domain  through  the  use  of  the  inverse  FFT. 
Frequency  domain  based  normalized  crosscorrelation  is 
described  by 
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p^(i)  =  It  u(^)-m):)  0.1.  N-i 

w-i  jv-i 

N  fc'O  N  k'O  (19) 

=  0  Otherwise 

where  X{k)  =  FFT  of  the  series  x{n) 

Y(k)  =  FFT  of  the  series  y  (n) 

*  =  conjugation 

Pj^(i)  =  normalized  correlation  function 
To  avoid  circular  correlation,  both  time  domain  pulse 
sequences  are  zero  padded  to  M,  data  points 

Ny  -  1 

where  =  number  of  data  points  in  x{n) 

Ny  =  number  of  data  points  in  y[n)  . 

To  take  advantage  of  the  processing  speed  of  FFT  algorithms, 
N,  is  made  a  power  of  two  by  increasing  its  length  with 
additional  zero  padding 


Ng  =  2‘"  Z  *  Ny  -  1 

where  m  is  an  integer  . 


Prior  to  any  zero  padding,  the  dc  component  is  removed  to 
create  zero  mean  pulse  sequences.  Both  sequences,  each  of 
length  N,  are  then  transformed  to  the  frequency  domain  through 
the  FFT  and  processed  using  the  frequency  domain  correlation 
technique.  A  MATLAB  implementation  of  frequency  domain 
crosscorrelation  is  given  in  i^pendix  B  (FreqCorr7 .m) . 
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C.  TINS  DirFBRXNCX  OF  ARRIVAL  (FOURIER  DOMAIN) 

The  normalized  croascorrelation  function  implemented  in 
the  frequency  domain  (Equation  19)  is  used  to  estimate  the 
TDOA.  The  croascorrelation  of  two  similar  signals  will  produce 
a  dominant  peak  in  the  output  of  the  correlator.  This  peak  is 
used  as  an  estimate  of  the  TDOA  between  the  two  signals.  For 
two  signals  that  have  not  been  corrupted  by  noise,  the 
location  in  time  of  this  peak  is  easily  determined.  Pulsed 
signals  that  have  been  distorted  by  noise  produce  a  noisy 
correlated  output.  The  output  noise  can  mask  the  TDOA  peak  in 
low  SNR  conditions.  Consequently,  the  probability  of  selecting 
a  noise  seunple  instead  of  the  true  TDOA  peak  increases  with 
decreasing  SNR. 

Figures  7,8  and  9  show  the  frequency  domain  based 
crosscorrelation  of  two  pulse  sequences  as  a  function  of 
decreasing  SNR.  A  reception  is  simulated  with  a  pulse  burst 
demodulated  by  one  receiver  fed  into  one  channel,  and  a  pulse 
burst  demodulated  by  a  second  receiver  fed  into  the  second 
channel  of  the  correlator.  Pulse  burst  one  has  a  transmission 
delay  of  zero.  Pulse  burst  two  has  a  delay  of  50  time  units. 
Therefore,  the  total  TDOA  is  50  time  units.  The  TDOA  peak  is 
clearly  seen  at  the  high  SNR  of  -i-7dB  and  +3dB.  At  OdB  SNR, 
noise  peaks  exceed  in  amplitude  the  true  TDOA  peak  causing  an 
incorrect  estimate  if  one  was  selected. 
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SNR  +7dB. 


SNR  +3dB. 

Visually  estimating  the  location  of  the  peak  with  the 
maximum  amplitude  in  the  output  of  the  correlator  is  a  easy 
task.  The  decision  rule  states  that  the  point  with  the  maximum 
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Figur*  9.  Normalized  frequency  domain  based  crosscorrelation. 
SNR  OdB. 


amplitude  is  selected  as  the  TDOA  estimate,  provided  the 
amplitude  exceeds  a  predetermined  threshold. 

When  many  blocks  of  data  must  be  processed,  visual 
inspection  of  the  correlation  output  is  not  possible.  Chapter 
V  discusses  the  design  of  a  threshold  for  automated  data 
processing. 
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IV.  NOISB  AMD  SIGNALS  IN  RBCBIVBR 


A.  NOISB  AMD  SIGNALS  IN  THB  RADAR  DBTBCTOR 

The  purpose  of  a  radar  receiver  is  to  process  and  detect 
the  reflected  energy  from  a  target  being  tracked  by  the  radar. 
An  I/Q  demodulator  is  assiomed  in  the  radar  receiver,  located 
after  the  last  intermediate  frequency  (IF)  amplifier  and  prior 
to  the  video  signal  display.  Data  to  be  processed  for  TDOA  can 
be  obtained  from  three  locations  on  the  I/Q  demodulator.  These 
locations  include  the  envelope  output,  the  envelope  sc[uared 
output  and  the  I/Q  channel  data. 

White  Gaussian  distributed  noise  is  assumed  to  be  added  to 
radar  pulses  during  their  tramsmission  through  the  atmosphere 
and  reception  by  the  receiver.  This  noise  contributes  to  the 
degradation  of  the  SNR  and  is  caused  by  fluctuations  in  the 
target  radar  cross-section  auid  the  loss  of  signal  energy  due 
to  the  propagation  distauice  between  the  target  auid  the 
receiver. 

At  typical  radar  frequencies  the  front  end  amplifier  of 
the  radar  receiver  contributes  to  this  noise  power.  The  noise 
power  for  a  thermally  noise  limited  receiver  is  defined  as  Nq, 
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No  =  k  B 

where  k  =  Boltzman's  constant 
B  =  IF  bandwidth 

(22) 

T  =  T  +T 
sys  «r  a 

=  receive  antenna  temperature 
r.  =  effective  noise  temperature 
of  the  receiver  . 


The  total  noise  is  modeled  as  zero  mean  Gaussian  distributed 
noise  with  a  probability  density  function  (pdf)  of 


Pjr(x) 


-  1  2o* 


(23) 


where  is  the  noise  variance  due  to  all  noise  sources  (i.e., 
receiver  front  end,  transmission  noise) . 

For  an  operational  radar,  either  a  pulsed  signal  plus 
Gaussian  noise  or  Gaussian  distributed  noise  alone  is  assumed 
at  the  input  of  the  I/Q  demodulator.  The  probability 
distribution  of  the  demodulated  signal  is  dependant  on  which 
of  the  three  outputs  of  the  I/Q  demodulator  is  used. 

1 .  Coherent  detection 

If  the  exact  frequency  eund  phase  of  the  received  pulse 
burst  is  known,  and  matches  the  frequency  and.  phase  of  the 
local  oscillator  in  the  receiver,  then  the  signal  can  be 
coherently  detected.  Coherent  detection  simplifies  the 
probeUsility  description  of  the  signal  amd  noise  at  the  output 
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of  the  I  and  Q  channels.  The  I  and  Q  sequences  have  an 
independent  Gaussian  distribution,  and  can  be  described  as 

1,0  -  N({nij{n)  ,mQ(n)}  ,a^)  (24) 

where  m^{n)  ,  irigin)  =  I  and  0  wean  values  respectively 

=  noise  variance  due  to  all  noise  sources  . 

Both  the  I  or  Q  channel  can  be  processed  using  the  correlation 
algorithm  followed  by  threshold  detection. 

Due  to  the  passive  nature  of  signal  reception  and 
demodulation,  the  exact  frequency  and  phase  of  the  received 
pulse  burst  are  not  available  at  the  receiver.  Therefore, 
coherent  detection  is  not  feasible  for  TDOA  estimation. 

2.  Envelop#  Squared  Detection 

For  a  zero  mean  Gaussian  distributed  noise  input  to 
the  I/Q  demodulator,  and  each  have  a  chi-squared 

distribution  with  one  degree  of  freedom.  Their  pdf's  can  be 
described  by 

Pr^y)  -  — j —  e  for  y  2  0 

a^2ny 

-  0  for  y  <  0  (25) 

where  o*  *  variance  of  I  and  0 
y  -  or  0^ 

The  envelope  squared  output  is  the  sum  of  and  and  is 
therefore  chi-squared  with  two  degrees  of  freedom.  This  pdf 
can  be  described  by 
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p^{z)  = 


for  z  i  0 


(26) 


e 

20^ 

where  z  =  +  0^ 

Envelope  squared  detection  is  attractive  because  the  Hilbert 
transform  can  be  used  to  enhance  the  TDOA  peak  in  the 
correlator's  output.  Unfortunately  in  a  low  SNR  environment 
the  noise  power  will  also  increase  as  the  noise  contribution 
is  squared. 

3.  Envelope  Detection 

A  common  form  of  receiver  demodulation  at  IF 
frequencies  is  the  envelope  detector.  For  all  simulations  in 
this  thesis /  envelope  detection  will  be  assumed.  The  output  of 
the  envelope  detector  is  the  modulation  envelope  of  the 
received  pulse  burst.  For  Gaussian  distributed  noise  at  the 
input  to  the  I/Q  demodulator,  the  noise  at  the  envelope  output 
is  Rayleigh  distributed  with  pdf  given  by 

Pu(u)  *  -^e  ,  u  >  0 

=  0,  otherwise  (27) 

where  *  variance  of  Gaussian  noise 
u  yjz  -  y/i^  *  0^  . 

For  a  pulsed  sinusoid  plus  Gaussian  noise  at  the  input  to  the 
I/Q  demodulator,  the  pdf  of  the  envelope  at  the  output  is 
Rician  distributed.  The  Rician  pdf  is  given  by 
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p«(r)  =  4®  ,  r  i  0 

=  0,  otherwise  (28) 

where  A  =  amplitude  of  received  sinusoid 
=  varia/2ce  of  Gaussian  noise 

where  lAz)  a  —  d0 

®  2n  Jo 

(29) 

and  r^Cz)  =  modified  Bessel  function  of  the 
first  kind  of  zero  order 

Figure  10  shows  a  block  diagraun  of  the  I/Q  demodulator  and  the 
various  pdfs  in  the  radar  receiver. 

B.  EMVBLOPB  COBBBIATIOM  OF  RAYLBIOT  MOISB 

Rayleigh  distributed  noise  is  produced  at  the  output  of 
the  envelope  detector  under  noise  only  conditions.  This  noise 
is  then  processed  by  the  correlator  during  the  TOOA  estimate. 
In  all  follow  on  discussions  and  simulations,  the  envelope 
processing  scheme  is  used. 

The  autocorrelation  of  Rayleigh  noise  as  discussed  earlier 
produces  a  large  triangular  bias  in  the  output  which  decreases 
as  the  time  lag  increases.  This  bias  is  produced  by  the  dc 
component  of  the  Rayleigh  distributed  noise. 

This  triangular  function  contains  no  information  relative 
to  the  TDOA  estimate.  The  large  numerical  values  in  the  bias 
can  potentially  limit  the  dyneunic  range  in  the  FFT  based 
correlation.  A  constant,  detection  threshold  cannot  be 
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I/Q  demodulator  in  the  radar  receiver. 


isg>lemented  at  the  output  of  the  correlator  with  this 
triangular  function  present.  For  these  three  reasons^  it  is 
advantageous  to  force  the  Rayleigh  or  Rician  distribution  to 
have  a  aero  mean.  This  will  remove  the  triangular  bias  formed 
by  the  correlation  process.  This  modification  of  the  received 
pulse  burst  does  not  degrade  the  TDOA  estimate.  For  all 
realizations  (i.e.,  blocks  of  data) r  the  dc  component  of  the 
Rayleigh  noise  will  be  removed. 
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1.  Shifted  Rayleigh  Hois* 


Forming  a  Raylalgh  distributed  time  series  x(n)  and 
subtracting  its  expected  value  produces  a  shifted  Rayleigh 
sequence.  This  algorithm  is  described  by 


Xi(n)  »  x(n) 


(30) 


where  in)  «  shifted  Rayleigh  noise  sequence 

o  *  Gaussian  noise  standard  deviation 


The  shifted  Rayleigh  time  series  has  a  zero  mezm.  The  new  time 
series  and  the  corresponding  time  domain  based  autocorrelation 
function  are  shown  in  Figure  11. 


Figure  11.  (a)  Shifted  Rayleigh  time  series,  (b)  Normalized 

shifted  Rayleigh  time  domain  based  autocorrelation. 


Histograsis  of  a  Rayleigh  pdf  and  a  shifted  Rayleigh 
pdf  are  shown  in  Figure  12.  Clearly ,  the  effect  of  removing 
the  dc  term  from  the  Rayleigh  distributed  time  series  creates 
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an  autocorrelated  output  that  does  not  have  the  limiting 
triangular  bias.  The  impulse  at  zero  time  lag  dominates  the 
correlation  output.  A  constant  detection  threshold,  at  least 
over  a  small  number  of  range  bins,  could  be  implemented  with 
this  t3^e  of  correlation  output. 

2 .  rrequeney  Oooain  Based  Correlation  of  Shifted  Rayleigh 
Noise 

The  frequency  domain  based  autocorrelation  of  x^Cn)  is 
now  exaunined.  Sequence  Xj^(n)  is  zero  padded  as  described  in 
Equation  21.  A  frequency  domain  series  is  created  using  the 
FFT  as  described  by  Equation  16.  The  autocorrelation  estimate 
is  derived  using  Equation  19  and  is  plotted  in  Figure  13.  The 
results  of  frequency  domain  based  autocorrelation  are  very 
similar  to  the  results  of  time  domain  based  autocorrelation 
using  shifted  Rayleigh  noise. 
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3.  Spectral  Xntaxpolation 

With  the  incorporation  of  the  FFT  in  modern  radar 
receivers,  it  is  possible  that  the  received  pulse  burst  has 
been  transformed  into  frequency  data  as  a  byproduct  of 
detection.  Processing  time  would  be  lost  in  converting  this 
data  back  into  a  time  series  to  be  correlated.  Of  course, 
there  is  no  guarantee  that  the  dc  component  has  been  removed. 
In  this  section  an  algorithm  is  given  to  perform  frequency 
domain  based  correlation  on  a  Rayleigh  distributed  noise 
sequence  while  avoiding  the  triangular  correlation  bias 
created  by  the  (possible)  dc  component. 

A  Rayleigh  distributed  noise  series  x(n),  N  points 
long,  is  transformed  into  an  ordered  sequence  of  N  complex 
coefficients  using  the  FFT  algorithm.  In  general,  a  sequence 
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X(k)  of  coefficients  is  indexed  from  zero  to  N-1  and  is 
described  by 

X{k)  =  FFT  {  x{n)  } ,  Oik  sN-1 
where  k  =  frequency  index  . 

The  FFT  creates  the  zero  indexed  term  X(0)  by  svimming  all  the 
elements  in  x(n).  X(0)  is  real  and  can  be  viewed  as  the  dc 
component  of  x<n).  By  removing  the  X(0)  coefficient  from  the 
frequency  sequence  we  have  in  effect  removed  the  dc  bias  from 
the  time  series  x (n) .  At  this  point,  the  modified  Fourier 
transform  X(k)  approximates  the  DFT  of  a  shifted  Rayleigh 
series . 

Time  domain  correlation  of  an  N  point  series  produces 
2N  data  points.  Frequency  domain  based  correlation  of  an  N 
point  sequence  produces  N  data  points.  To  force  frequency 
domain  based  correlation  to  equal  time  domain  based 
correlation  (i.e.,  avoid  circular  correlation),  the  number  of 
terms  in  the  frequency  series  must  be  doubled.  The  frequency 
series  is  expanded  by  storing  each  pair  of  complex  frequency 
coefficients  at  a  location  twice  the  value  of  the  original 
index . 

As  a  first  approximation,  the  data  points  between 
adjacent  coefficients  are  linearly  interpolated  in  the 
expanded  array.  Figure  14  illustrates  the  effect  of  FFT 
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interpolation  and  the  zeroing  of  the  dc  terns.  Correlation 
results  as  described  by  Equation  19  are  shown  in  Figure  15. 


Re(  X(Kl  }  or  im{  Xlkj  f 

A 

x(r,i/2) 

N/2 

I  =  inte'ooiated  coefficient 
X(i>i  =  Fojrier  coefficient 


Figure  14.  FFT  interpolation,  and  zeroing  of  the  dc 
coefficient . 
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Figure  15.  Normalized  frequency  autocorrelation  through 
interpolation . 
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Clearly,  the  triangular  bias  in  the  correlation  output  has 
been  removed.  A  small  impulse  at  a  time  lag  of  -N  is  observed, 
which  can  be  disregarded. 

For  block  processing,  typically  only  correlation 
products  within  ±10%  of  the  zero  time  lag  position  are 
considered  to  be  statistically  reliaOsle.  Therefore,  the  TDOA 
estimate  is  formed  within  this  correlation  band.  The  impulse 
at  -N  falls  outside  of  the  TDOA  evaluation  area  and  will  not 
cause  a  false  detection. 

C.  CORRELATION  OF  PULSES  IN  ADDITIVE  RAYLEIGH  NOISE 

Time  domain  based  correlation  of  a  pulse  train  in  additive 
Rayleigh  noise  should  produce  the  same  result  as  frequency 
domain  based  interpolation  followed  by  am  equivalent 
correlation.  A  comparison  was  performed  using  a  50%  duty  cycle 
pulse  train  in  additive  shifted  Rayleigh  noise. 

Time  correlation  was  performed  on  the  pulse  burst.  The  dc 
component  of  the  noisy  pulse  train  is  removed  forming  a  zero 
mean  signal.  The  time  domain  based  autocorrelation  of  this 
sequence  is  shown  in  Figure  16(a). 

The  pulse  burst  was  also  transformed  into  the  frequency 
domain  using  the  FFT  with  the  dc  component  intact.  N  Fourier 
coefficients  are  expanded  and  the  X(0),  X(±l)  coefficients 
zeroed  out.  The  sequence  is  correlated  after  interpolation 
using  Equation  19  and  plotted  in  Figure  16(b). 
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autocorrelation.  Zero  time  lag  between  each  channel. 

Clearly,  the  two  correlation  algorithms  produce  similar 
outputs.  Both  have  a  dominant  pe-»lc  in  their  output  for  TDOA 
estimation.  Frequency  domain  based  correlation  displays  a 
exponential  decay  of  the  sidelobes  about  the  main  peak.  This 
effect  is  not  observed  in  time  autocorrelation,  and  is  caused 
by  the  first  order  interpolation  used  in  the  frequency  domain 
based  method.  A  more  robust  interpolation  algorithm  will 
eliminate  the  exponential  decay,  yielding  the  desired  linear 
decay. 

Next,  time  and  frequency  domain  based  crosscorrelation  are 
performed  on  two  pulse  sequences.  One  pulse  sequence  lags  the 
other  by  20  units.  Both  time  and  frequency  domain  based 
crosscorrelation  of  the  two  sequences  are  plotted  in  Figure 
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17.  Clearly,  both  correlation  algorithms  produce  a  peak  at  20 
units.  The  nonlinear  decay  of  the  sidelobes  surrounding  the 
dominant  peak  is  again  seen  in  frequency  domain  correlation. 


Figure  17.  (a)  Normalized  time  domain  based  crosscorrelation, 
(b)  Normalized  frequency  domain  based  interpolation  and 
crosscorrelation.  20  unit  time  lag  between  each  channel. 


In  summary,  the  frequency  domain  based  interpolation  and 
crosscorrelation  algorithm  produces  an  output  comparable  to 
time  domain  based  crosscorrelation.  It  offers  a  method  to 
correlate  pulse  bursts  collected  in  the  frequency  domain. 
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V.  THRESHOLD  DESIQI  FOR  EHVBLOPB  CORRELATION 


A.  INTRODUCTION 

A  threshold  at  the  output  of  a  correlator  is  used  to 
estimate  the  location  of  a  TDOA  peak  in  a  constant  false  alarm 
rate  (CFAR)  receiver.  The  value  of  the  threshold  is  a  function 
of  the  correlation  output  mean  and  variance. 

B.  CORRELATION  OUTPUT  AND  CONSTANT  VARIANCE 

If  the  variance  of  the  correlator,  when  driven  by  two 
signals  at  the  input,  maintains  a  fixed  value  2d>out  some  mean, 
a  constant  threshold  above  the  meam  can  be  used.  This  allows 
the  design  of  a  constant  false  alarm  rate  detector.  The 
largest  correlation  peak  that  also  crosses  the  threshold  is 
defined  as  the  TDOA  estimate.  The  constant  threshold  design  is 
the  easiest  to  in^lemeac. 

C.  CORRELATION  OUTPUT  AND  TRIANGULARLY  SHAPED  VARIANCE 

The  crosscorrelation  of  two  zero  mean  noise  sequences  with 
the  same  number  of  data  points,  produces  an  output  whose 
variance  has  triangular  form.  Hence,  a  constant  detection 
threshold  is  not  easily  implemented.  Instead,  a  threshold  must 
be  designed  that  has  a  slope  that  matches  the  slope  of  the 
correlation  output.  The  design  of  a  sloping  threshold  is 
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difficult,  particularly  if  the  slope  of  the  correlation  peak 
changes  with  time  (i.e.,  changing  noise  environment). 

A  solution  is  to  mathematically  force  the  triangular 
correlation  variance  to  be  flat.  This  is  done  by  using  a 
bowtie  correction  of  the  form  1/(N-|(!  |  A  bowtie  correction 
normalizes  the  correlation  output  to  have  a  constant  (lag 
independent)  variance.  This  design  suffers  from  time  changes 
in  the  correlation  output  forcing  the  bowtie  correction  to  be 
recalculated.  Note  that  the  aunplitude  of  the  TDOA  peak  will 
itself  be  dependant  on  the  type  of  correction  to  the 
correlation  output. 

D.  TBBKSHOLD  BASBO  ON  ZIRO  LAG 

The  zero  lag  threshold  method  uses  the  variance  of  the 
zero  lag  correlation  output  to  calculate  a  constant  detection 
threshold.  This  is  indicated  in  Figure  18. 

The  correlation  output  becomes  statistically  less  reliable 
the  further  removed  the  correlation  products  are  from  the  zero 
time  lag  position.  To  use  correlation  intelligently,  the 
correlation  output  is  typically  examined  ±10%  of  the  distance 
from  the  zero  time  lag  position.  Consequently,  only 
correlation  output  ±0.1N  from  zero  time  lag  will  be  used  in 
TDOA  estimate  simulations. 

The  detection  method  in  this  thesis  uses  a  constant 
threshold  against  the  correlation  output.  The  fixed  threshold 
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was  selected  because  the  correlation  output  variance  is,  to  a 
first  approximation,  nominally  flat  within  the  tO.lN  boxindary. 

The  detection  scenario  that  drives  this  thesis  defines  a 
minimum  radar  pulse  width  of  one  microsecond  (p.sec)  and  a 
burst  length  of  four  millisecond  (msec)  .  A  50%  duty  cycle 
radar  pulse  is  assumed.  There  are  2000  pulses  per  burst,  which 
when  correlated,  will  generate  an  output  that  is  8000)isec 
(±4000fisec)  in  length.  If  ±10%  of  the  correlation  output 
around  zero  time  lag  are  compared  to  the  threshold  value, 
±400p.sec  will  be  tested.  For  radar  pulses  traveling  at  the 
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speed  of  light,  a  maximum  distance  of  ±120]cm  is  therefore 
examined  to  form  a  TDOA  estimate. 

Correlation  values  at  lags  -O.IN  to  0 . IN  will  be  compared 
in  magnitude  to  the  zero  time  lag  based  calculated  threshold 
value . 


B.  THRESHOLD  CALCDIATION  FOR  CFAR 

The  detection  threshold  is  calculated  using  the  output 
terms  of  the  crosscorrelation  algorithm  when  noise  only  is 
injected  into  the  correlator.  For  convenience.  Equation  8  is 
repeated  here, 

w-UI-i 

=  E  yi*i 

(32) 

and  are  zero  mean  Rayleigh  (shifted) 
distributed  noise  sequences  . 

The  correlation  function  at  all  lags  of  interest  is  thought  to 
be  Gaussian  distributed.  According  to  Equation  32,  every 
estimate  is  the  sum  of  a  fairly  large  number  of  products 
(i.e.,  central  limit  theorem).  The  average  value  of  the 
correlation  output  between  the  ±0.1N  endpoints  is 


—  E 


r^(i) 


where  a  =  number  of  data  points  between  the 
0  and  0 .  IN  lag  posi  tions  in  the 
correlation  output  . 


(33) 


The  variance  of  r^^d)  is  calculated  from 
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Ot  = 


2a  +  1 


i*-a 


(r,y{l)  - 


E  ) 


(34) 


For  Gaussian  distributed  noise  the  probability  of  false  alarm 
(Pf.)  is  a  function  of  the  detection  threshold 


Pfa 


where  T  =  detection  threshold 


(35) 


For  a  given  P,.,  the  threshold  can  be  calculated  from 

T  =  erfc'^  {PfJ  Ot 


(36) 

where  erfc  =  complementary  error  function  . 

Equation  36  relates  the  threshold  to  the  Gaussian  distributed 
noise  variance. 


F.  VERIFICATION  OF  CFAR  TBRISROLD  AND  ROC  CDRVX 

Using  the  above  equations,  thresholds  were  calculated  for 
the  time  correlation  detector  given  in  Equation  32  for  three 
different  Pfa's.  A  MATLAB  simulation  was  performed  to  compare 
the  measured  false  alarm  rate  of  the  time  correlation  detector 
to  the  theoretical  false  alarm  rate. 

1 .  F^a  Simulation 

Two  100  point  sequences  of  zero  mean  Rayleigh  noise 
were  injected  into  the  correlation  detector.  The  output  of  the 
correlation  detector  between  the  -O.IN  and  -i-O.lN  lag  points 
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was  compared  to  the  calculated  threshold  to  measure  the  false 
alarm  rate.  The  number  of  points  in  the  correlation  detector 
output  that  exceeded  the  threshold  for  a  given  P,.  were 
recorded.  The  experimental  P,,  is  given  as  the  ratio  of 
improper  threshold  crossings  to  the  total  ntamber  of  monitored 
range  cells  of  the  correlator  output  (see  Table  1) .  For  each 
P,.,  500  realizations  were  performed.  Both  theoretical  and 
experimental  results  are  listed. 


TABLE  1 

THEORETICAL  PROBABILITY  OF  FALSE  ALARM 

Simulation  # 

P„=0.001 

1 

0.0100 

0.00076 

0 

2 

0.0099 

0.00090 

0 

3 

0.0110 

0.00085 

• 

The  data  in  Table  1  show  that  a  detection  threshold 
can  be  established  for  a  correlation  detector,  using  the 
assumption  of  a  Gaussian  distributed  output  for  Pf.' s  of  0.01 
and  0.001.  For  the  lower  P,.  value  of  0.0001,  the  threshold 
calculation  fails,  but  in  a  positive  sense.  The  measured  P,. 
of  0  is  below  the  theoretical  Pf.  of  0.0001.  Either  the 
correlation  variance  is  smaller  than  Equation  36  predicts,  or 
more  terms  must  be  suxnmed  in  the  correlation  output  to  better 
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approximate  the  Gaussian  pdf  as  indicated  by  the  central  limit 
theorem. 

2 .  TDOA.  Siotulation 

The  output  of  the  envelope  detector  will  either  be 
Rayleigh  or  Rician  (non-central  Rayleigh)  distributed.  In 
those  portions  of  the  data  where  there  is  no  signal  energy, 
the  output  will  be  Rayleigh  distributed.  Where  the  signal  is 
present,  the  envelope  detector  output  will  be  Rician 
distributed . 

The  purpose  of  the  time  domain  based  correlation 
detector  is  to  measure  the  TDOA  of  pulsed  sequences  embedded 
in  Rayleigh  noise.  Noise  in  low  SNR  environments  will  mask  the 
TDOA  peak.  A  MATLAB  simulation  was  performed  to  exaunine  the 
performance  of  the  correlation  detector  in  a  decreasing  SNR 
environment . 

Two  pulse  bursts  embedded  in  additive  Rayleigh  noise 
were  created.  The  second  pulse  burst  lagged  the  first  by  three 
time  units.  Each  pulse  burst  is  a  sequence  of  ten  pulses,  each 
having  a  period  of  ten  units.  Each  pulse  in  the  burst  had  a 
50%  duty  cycle. 

The  SNR  established  the  amplitude  relationship  between 
the  noise  luid  signal  sequences.  First,  two  sequences  of  zero 
mean  unit  variance  Gaussian  distributed  noise  were  created. 
Rayleigh  noise  was  generated  from  the  two  Gaussian  sequences. 
The  Rayleigh  noise  had  a  mean  of  1.2533  and  a  variance  of 
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0.4292.  The  theoretical  second  moment  of  the  Rayleigh  noise  is 
two  and  was  verified  in  the  MATLAB  simulation.  The  second 
moment  is  defined  as  the  average  power  of  the  Rayleigh  noise. 
The  amplitude  of  the  nonzero  elements  in  each  pulse  burst  was 
then  scaled  to  meet  the  desired  SNR,  using  the  measured 
average  power  of  the  Rayleigh  noise. 

A  detection  threshold  was  calculated  based  on  shifted 
Rayleigh  noise  using  Equation  36  where  was  assumed  to  be 
known  (i.e.,  simulation  information).  In  practical 
applications,  would  be  estimated  from  the  data.  The  SNR  of 
the  two  bursts  were  varied  from  OdB  to  18dB  in  the  different 
simulations.  For  each  SNR,  the  two  sequences  were  correlated 
using  Equation  32. 

Graphs  of  the  correlation  output  for  a  SNR  of  IdB  are 
given  in  Appendix  C  for  Pf,'s  of  0.01,  0.001,  0.0001  and 
0.00001.  The  threshold,  which  is  a  function  of  the  Pf.  and 
noise  vari2mce,  is  shown  to  increase  with  decreasing  P,.. 

The  correlation  peak  at  a  time  lag  of  three  was  the 
correct  TDOA  estimate.  If  the  maxlmvun  correlation  peak  in  the 
output  was  the  third  lag  point,  a  correct  TDOA  estimate  was 
obtained.  For  each  SNR,  30  TDOA  estimates  were  made.  The 
correlation  fxinction  output  100  data  (lag)  points  per  TDOA 
estimate.  The  TDOA  estimate  was  made  on  ±10  data  points  on 
each  side  of  the  zero  time  lag  position  in  the  correlation 
output.  The  percent  correct  TDOA  estimate  for  each  SNR  was 
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calculated  and  is  presented  as  a  receiver  operating 
characteristic  (ROC)  curve  in  Figure  19. 

The  ROC  curve  shows  that  as  the  SNR  increases,  the 
performance  of  the  correlation  detector  increases.  The 
performance  of  the  correlation  detector  in  Figure  19  can  be 
improved  if  multiple  sequential  looks  at  the  same  SNR  were 
allowed. 
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00 


Figur*  19.  Percentage  of  estimated  TDOA  correct  versus 
SNR  for  time  domain  based  correlation  detection. 
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VI.  FFT  BRROR  ANALYSIS 


A .  INTRODUCTION 

This  chapter  discusses  several  sources  of  processing 
errors  occurring  within  the  FFT,  and  the  effect  of  these 
errors  on  the  FFT  output  signal-to-noise  ratio. 

The  computation  of  a  FFT  using  fixed  point  arithmetic, 
generates  internal  processing  errors.  These  errors  are  due  to 
the  fixed  niimber  of  bits  in  the  trigonometric  look  up  table, 
the  truncation  of  arithmetic  results  during  addition  and  the 
scaling  of  results  during  the  individual  butterfly 
calculations . 

The  error  sources  can  be  modeled  as  independent  additive 
noise  components  at  the  output  of  the  FFT.  Their  net  effect  is 
to  reduce  the  ideal  SNR  at  the  output  of  the  FFT. 

Each  of  these  processing  error  sources  will  be  described 
below.  The  material  in  this  chapter  borrows  from  [Ref.  7] . 

B.  COSIMB/SIMB  TABLB  NOISB 

The  computation  of  a  butterfly  algorithm  in  the  FFT 
involves  multiplication  by  the  complex  coefficient  (twiddle 
factor) 
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(37) 


e  =  coai^m)  -  j  sin(-^/n) 
N  N 

=  W”  . 

A  FFT  butterfly  is  shown  in  Figure  20. 


Figure  20.  FFT  butterfly  calculation. 

These  trigonometric  coefficients  are  precalculated,  quantized 
to  (B)  bits  (where  B  is  the  size  of  the  con^uter  word 
including  sign  bit) ,  and  stored  in  memory  for  future  table 
loo)cup  during  the  FFT  processing.  The  coefficients  N**  ■  1  and 
a  -  j  can  be  obtained  in  an  exact  form  and  therefore  do  not 
contribute  to  the  quantization  noise.  Quantization  of  the 
other  (twiddle)  coefficients  stored  in  the  sine  and  cosine 


50 


taUsle  introduce  noise  into  the  output  of  the  FFT.  The  noise 
power  created  by  truncation  of  the  sine  and  cosine 
coefficients  is  defined  by  [Ref.  7]  to  be 


(38) 

where  =  trigonometric  noise  power 

B  =  number  bits  in  computer  word  . 


C .  TRUMCHTIOM  NOXSB 

The  DFT  of  a  finite  duration  sequence  {x(n) }  is  defined  in 
Equation  18  and  is  repeated  here  for  convenience 

-i  ^*ji 

X{k)  =  x{n)  e  ^  ,  0  i  ic  i  N~l  (39) 

n«0 

The  product  formed  requires  four  real  multiplications  because 
both  the  exponential  and  the  input  data  are  complex  numbers, 
B  bits  in  length.  Each  multiplication  can  produce  a  result  of 
2B  bits  which  must  be  truncated  from  2B  to  B  bits,  hence  there 
are  four  quauntization  errors  for  each  complex  valued 
multiplication.  The  four  quantization  errors  can  be  described 
as  four  independent  noise  sources.  The  truncation  noise  power 
is  defined  by  [Ref.  7]  as 
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(40) 


where  =  truncation  noise  power 

B  =  number  of  bits  in  product  after  truncation 


D .  SCALING  MOZSB 

The  FFT  as  an  algorithm  processes  a  data  vector  of  length 
N  by  successive  passes  at  the  vector.  During  each  pass,  the 
algorithm  performs  N/2  butterfly  operations.  Each  butterfly 
retrieves  two  complex  niimbers  from  memory,  performs  the 
butterfly  computation  and  returns  two  complex  numbers  to  the 
saune  memory  address.  The  complex  numbers  returned  to  memory  by 
the  butterfly  can  have  a  greater  magnitude  than  the  two 
complex  numbers  initially  fetched  from  memory  by  the 
butterfly.  In  order  for  the  results  of  the  butterfly  operation 
to  fit  in  the  fixed  word  length  of  the  memory  of  the  computer, 
the  results  must  be  truncated  (scaled) .  Scaling  is  performed 
by  dividing  by  two  the  entire  data  array,  and  is  implemented 
in  software  by  a  right  shift  and  an  increment  of  the  arrays 
exponent  register. 

During  the  (m+1)  pass  of  the  data,  the  butterfly  algorithm 
selects  two  data  points  A(m)  emd  B(m)  and  returns  to  memory 
A(m-<-l)  auid  B(m'fl).  The  truncation  error  results  from  the 
addition  of  two  numbers  of  like  sign  in  the  upper  part  of  the 
butterfly  algorithm  or  the  stjbtraction  of  two  numbers  of 
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opposite  sign  in  the  lower  part  of  the  butterfly  algorithm. 
The  input  to  the  butterfly  is  fixed  to  have  maximtjm  real  and 
imaginary  values  of  ±1 .  The  output  of  the  butterfly  is  fixed 
to  have  a  maximum  value  of  ±2,  double  the  input  value.  To 
prevent  overflow  of  the  butterfly  output,  either  a  prescale  by 
1/2  is  required  prior  to  entering  the  butterfly  or  an  extra 
bit  must  be  availzUsle  in  memory  to  accommodate  the  possible 
gain  of  two. 

One  method,  automatic  prescaling,  truncates  the  least 
significant  bit  of  the  butterfly  output  and  represents 
processing  noise.  In  the  case  for  which  no  scaling  is 
necessary  for  the  pass,  this  truncation  represents  significant 
processing  error  (noise) . 

A  second  method,  data  dependent  scaling,  is  performed  only 
if  a  butterfly  in  the  data  pass  overflows  without  scaling.  If 
any  word  in  memory  has  an  overflow  bit  set,  all  words  are 
shifted  right  as  they  are  fetched  from  memory  for  the  next 
pass  of  the  FFT.  The  total  number  of  right  shifts,  p,  executed 
during  the  FFT  process  is  used  at  the  output  of  the  FFT  as  a 
scale  factor  of  2’’.  This  scale  factor  is  used  at  the  output  of 
the  FFT  so  that  the  total  processing  gain  is  maintained. 
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The  noise  power  created  by  the  scaling  of  the  data  is 


defined  by  [Ref.  7] 


^23  4 


wheze  o|  -  scaling  noise  power 


(41) 


B.  OUTPUT  SIGNXL  TO  MOISB  RXTXO 

Figure  21  shows  an  error  model  of  a  butterfly  in  an  FFT. 


q^(fri) 


Figure  21.  FFT  error  model. 
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The  input  data  terms  to  a  butterfly  can  be  described  by 


A{m)  =  Aim)  +  Nj^im) 

Sim)  -  Bim)  +  Ngim) 

(42) 

where  Nj^,  =  additive  noise  terms 
Aim)  ,  Bim)  =  input  data  terms 
Aim)  ,  Sim)  =  noise  corrupted  data  . 

Without  the  effects  of  noise  the  butterfly  algorithm  produces 

two  outputs  described  by 

A(fl7+l)  =  Aim)  +  wim)Bim) 

Bim*l)  =  Aim)  -  Wim)Bim)  (43) 

where  Wim)  =  twiddle  factor  . 

The  error  due  to  the  A  component  will  be  calculated.  The  error 
due  to  the  B  component  is  statistically  equivalent  to  the 
error  in  the  A  component.  Including  the  effects  of  noise,  the 
output  of  the  butterfly  cam  be  described  by 

A(in+l)  *  [Aim)  ^^J^im)  +  [Bim)  *  Ngim)  ^q^im)]- 

[W(m)  +  <3r»,(m)]  +  Qjim) 

where  g^im)  =  scale  noise 

q^im)  =  sine/ cosine  table  noise 
q^im)  =  truncation  noise  iin  multiply) 

(44) 

The  total  power  at  the  output  of  the  butterfly  is  then 
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(45) 


E{A{m*l)  A*(/n+l)}  =  E[A(m)A*{m)}  E{B(.m)  B'  iw) }  + 

E{Nj^im)NA(m)]  *  E{NB(in)N;{m)]  + 
ElQsim)  gs(.m) }  +  (m) }  + 

E{3{m)  B^  {m)  ]  E{g^(m)  gC(m)  }  + 
E{gj.(m)  gT(Tn)  }  . 


Redefining  the  terms  of  the  total  power  leads  to 


E{A{m*l)A'  (m+1) } 

where  {m) 
olim) 
a%{m) 
olim) 


im)  +  (m)  +  Na  (m)  +  Ng  (m)  + 

2a\{m)  +  S^{m)a%{m)  +  o|(/n) 

input  signal  power 
scaling  noise  power 
truncation  noise  power 
sine/ cosine  table  noise  power 


Combining  like  terms  results  in 


£{^(m+l)^*  (m+1) }  =  2S^{m)  *  2Nl[m)  +  2al{m)  + 

S^{m)al{m)  *  olim)  . 


(47) 


This  equation  describes  the  total  power  as  a  function  of  the 
signal  power  and  the  individual  independent  uncorrelated  noise 
terms.  The  total  power  output  from  the  butterfly  can  also  be 
described  as  a  fxinctlon  of  its  input  signal  power  eind  the 
total  additive  noise  power  generated  within  the  butterfly  as 

£{^(in+l)i[* (xn+l) }  =  £{A(in+l) A* (m+l) }  +  #4g» 

Redefining  these  terms 

E{Aim*l)A*  im*l)}  *s2(m+l)  (49) 
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Equation  47  and  49  are  equivalent.  The  signal  and  noise  power 
for  the  (m+1)  pass  can  be  formed  from  those  terms  in  the 
proceeding  pass 

SM/n+1)  +  =  2S^(m)  +  2Nl{m)  +  2a\[m)  + 

(m)  a%(m)  +  o%{m)  . 

Equating  the  signal  terms  in  Equation  50  yields 

S2 (m+l)  =  2S^  (m)  . 

Equating  the  noise  terms  in  Equation  50  yields 

=  2Nl(m)  +  2a\(m)  +  3^  (m)  o]r(m)  +  a%(m)  . 

Equation  51  can  be  rewritten  to  show  how  the  input  signal 
power  increases  as  a  function  of  m  as  it  passes  through  the  m^'’ 
stage  of  a  butterfly  of  an  FFT 

(53) 

here  we  modeled  5^(0)  =  S^. 

The  signal  power  input  to  the  butterfly  is  defined  as 

52  =  .  (54) 

9 

Rewriting  Equation  52  to  include  this  exponential  increase  of 
S*  (m)  gives 


(50) 


(51) 


(52) 
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=  2Nl(m)  ^  2al{m)  +  2‘^S^olim)  +  Or^m)  .  (55) 

Equation  55  shows  the  relationship  between  the  noise  power  of 
a  butterfly  stage  as  a  function  of  the  noise  power  from  the 
proceeding  butterfly  stage.  The  noise  and  signal  power  can  now 
be  calculated  using  Equations  53  and  55  for  m  passes  through 
a  butterfly.  Once  calculated,  the  signal  to  noise  ratio  can  be 
formed. 

F.  RIGHT  JUSTIFIED  DATA 

The  length  of  a  word  in  the  computer  is  defined  to  be  16 
bits  (B)  and  accepts  a  12  bit  (b)  data  word.  Right  justified 
data  places  the  least  significant  bit  (LSB)  of  the  b  bit  data 
word  in  the  LSB  of  the  B  bit  processor  word.  The  output  noise 
to  signal  ratio  is  now  calculated  for  ten  passes  through  a 
butterfly  using  right  justified  data.  The  following 
assumptions  are  used  in  the  calculation; 

1.  The  input  signal  is  incoherent. 

2.  Scaling  is  performed  every  other  pass  after  the  most 
significemt  bit  (MSB)  of  data  reaches  the  next  to  MSB 
of  the  processor. 

3.  The  first  two  passes  do  not  use  the  trigonometric 
table. 

After  the  tenth  pass  (i.e.,  FFT  size  of  1024),  it  is  shown 
(i^pendix  D)  that  the  noise  to  signal  power  is 
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(56) 


Aal  * 

4  51  52 


B  and  b  are  previously  assumed  bo  be  16  and  12  bits  in  length 
respectively.  Substituting  these  values  into  Equation  56 


_7  _-32  1  o  -30 

1  12^^  +  1  24 

4  J.2-8  51  J,2-e 

9  9 


.  (57) 


The  noise  to  signal  ratio  reduces  to 


-67.5  db  -95.1dB  -71.1  dB  -79.1  dB  (58) 

input  trig  truncation  scaling 

noise  noise  noise  noise 


Clearly,  the  dominant  degradation  noise  is  the  truncation 
table  noise  (-71.1dB),  followed  by  the  scaling  noise  to  signal 
ratio  (-79.1dB)  .  For  right  justified  data,  the  signal  to  noise 
ratio  is  calculated  to  be  65.7ldB. 


G.  van  JUSTXFZSD  xaax 

Left  justified  data  places  the  most  significant  bit  (MSB) 
of  the  b  bit  data  word  in  the  MSB  of  the  B  bit  FFT  processor 
word.  The  output  noise  to  signal  ratio  is  now  calculated  for 
10  passes  through  a  butterfly  for  left  justified  data.  The 
following  assumptions  are  made 
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1.  The  signal  is  incoherent. 

2.  Data  is  left  justified  on  input  and  scaling  is 
performed  every  other  pass  starting  with  the  first 
pass . 

3 .  The  first  two  data  passes  do  not  use  the 
trigonometric  table. 

It  cam  be  shown  that  after  the  tenth  pass  through  the 
butterfly  that  the  noise  to  signal  ratio  is 


^  1 


4  52 


(59) 


B  and  b  were  previously  assumed  to  be  16  and  12  bits 
respectively.  Substituting  these  values  into  Equation  59  we 
obtain 


1.  ^  0-32 

4  12 


(60) 


The  noise  to  signal  ratio  reduces  to 


^(10) 

52 


-61 .5  dB  -95.1  dB  -89.1  dB  -77 . 1  dB  (61) 

input  trig  truncation  scaling 

noise  noise  noise  noise  . 


Clearly,  the  worst  noise  contribution  to  the  SNR  is  the 
scaling  noise  (-77.1dB),  followed  by  the  truncation  noise  (- 
89.1dB) .  For  left  justified  data,  the  signal  to  noise  ratio  is 
calculated  to  be  67.0dB. 
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Further  analysis  shows  that  for  large  FFT's  (i.e.,  in  the 
range  of  2“  and  above) ,  the  trigonometric  table  noise  is  the 
dominant  noise  source.  A  simulation  was  performed  using  an 
AMDAHL  machine  to  examine  the  error  generated  in  performing  a 
FFT  and  an  iFFT.  Figure  22  shows  the  experimental  setup. 


AMDAHL  Machine 


Figure  22.  FFT  processing  error  flowchart  [Ref.  8] . 

The  error  that  was  generated  was  tabulated.  Figure  23 
illustrates  the  error  as  a  function  of  FFT  length. 

B.  SUMMABY 

For  either  left  or  right  justified  data,  the  domineint 
noise  at  the  output  of  the  1024  point  FFT  is  the  scaling  or 
truncation  noise.  The  SNR  at  the  output  of  the  FFT,  for  left 
justified  data  is  €7.0dB  and  exceeds  the  SNR  for  right 
justified  data  which  is  65.71dB.  Clearly,  left  justified  data 
maximizes  the  con^utational  SNR  at  the  output  of  the  FFT  for 
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MAXIMUM  ERROR  IN  SINGLE 
PRECISION  FFT 


FFT  LENGTH 


Figur*  23.  FFT  error  as  a  fiinction  of  transform  length 


[Ref.  8]. 
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the  values  of  b,  B  and  N  used  in  the  above  excunples . 


Comparing  Equation  56  to  Equation  63,  we  see  that  by 
increasing  the  transform  size  from  to  2^^  increases  the 
noise  to  signal  ratio.  The  trigonometric  table  noise  is 
increasing  faster  than  either  the  scaling  or  truncation  noise. 
The  signal  to  noise  ratio  is  reduced,  from  65.71dB  to  65.68dB. 
As  the  trams  form  size  is  increased  to  accommodate  larger  data 
sets,  the  trigonometric  table  noise  becomes  the  dominant  error 


term. 


VTI.  SUMMARY 


A.  CONCLUSIONS 

A  frequency  domain  based  algorithm  was  developed  and 
tested  to  estimate  the  differential  arrival  time  of  a  pulsed 
radar  signal  collected  by  two  passive  sensors.  For 
convenience,  the  actual  implementation  was  performed  through 
time  domain  processing  even  though  frequency  domain  processing 
is  advocated.  The  performance  of  the  algorithm  is 
characterized  by  a  ROC  curve  as  a  function  of  SNR.  For  3dB  SNR 
and  a  P,,  of  0.01,  the  probability  of  making  a  correct  TDOA 
estimate  exceeds  that  of  an  incorrect  estimate  for  a  single 
look.  At  18dB  SNR,  the  probability  of  making  a  correct  TDOA 
estimate  is  100  percent  for  all  Ffa's.  If  multiple  looks  (i.e., 
multiple  bursts)  are  allowed,  the  probability  of  making  a 
correct  TDOA  decision  at  each  SNR  will  increase. 

An  I/Q  demodulator  is  assumed  in  the  radar  receiver.  The 
pdf  of  the  signal  driving  the  correlation  detector  is 
determined  from  where  it  is  obtained  in  the  I/Q  receiver.  The 
pdf,  depending  on  selection,  will  be  zero  mean  Gaussian  versus 
non-zero  mean  Gaussian,  or  chi-squared  versus  non-central  chi- 
squared,  or  Rayleigh  versus  Rician  (non-central  Rayleigh) . 
This  thesis  assumes  an  envelope  detector  at  the  output  of  the 
I/Q  demodulator.  The  envelope  detector  output  has  a  Rayleigh 
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or  Rician  pdf  and  drives  the  correlation  detector.  The 
calculation  of  a  CFAR  threshold  assumes  a  Gaussian  pdf  at  the 
output  of  the  correlation  detector.  For  small  time  lags, 
summation  of  the  output  of  the  correlator  may  produce  a 
probability  distribution  that  deviates  from  Gaussian.  This 
deviation  would  produce  a  biased  threshold  calculation.  To 
minimize  the  bias,  a  sufficient  number  of  terms  must  be  summed 
at  the  output  of  the  correlator  to  allow  the  Gaussian 
approximation . 

If  the  received  pulses  are  collected  in  the  frequency 
domain,  a  spectral  domain  based  correlation  algorithm  can  be 
implemented.  An  algorithm  is  given  in  this  thesis. 

For  any  digital  signal  processing  algorithm  that  uses  the 
FFT,  processing  errors  must  be  considered.  For  large  transform 
sizes  the  trigonometric  noise  power  dominates.  The  length  of 
the  trigonometric  coefficient  word  affects  the  degradation  of 
the  SNR  at  the  output  of  the  FFT.  The  larger  this  word  the 
smaller  the  noise  power.  For  a  given  transform  size,  left 
justified  data  will  have  a  higher  output  SNR  than  right 
justified  data. 

For  a  correct  TDOA  estimate,  the  correlation  algorithm 
requires  the  reception  of  the  leading  pulses  in  the  radar 
pulse  burst.  If  the  pulses  are  received  in  an  adequate  SNR, 
the  correlation  algorithm  will  produce  a  well  defined  peak 
allowing  the  TDOA  estimate.  Should  the  environment  degrade  the 
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destructive  interference  due  to 


received  signal  (i.e., 
multipathing,  weather  or  terrain) ,  the  correlation  algorithm 
should  be  discarded  in  favor  of  the  traditional  angle  of 
arrival  (AOA)  algorithm  (i.e.,  maximize  received  energy). 

If  the  radar  uses  a  staggered  PRF,  the  reception  of  the 
first  pulse  or  several  pulses  has  little  impact  on  the 
position  of  the  correlation  peak.  Under  this  condition,  the 
TDOA  algorithm  is  superior  to  the  AOA  algorithm. 

B.  RBCOMMBMDATIONS 

A  ROC  curve  should  be  constructed  for  the  correlation 
algorithm  used  in  this  thesis  for  multiple  TDOA  looks.  An 
intensity  display  could  be  designed  that  would  plot  as  a 
fxinction  of  time  (i.e.,  snapshots),  those  lag  points  chosen  as 
TDOA  estimates.  Patterns  displayed  on  the  intensity  plots 
would  allow  the  user  to  visually  determine  the  correct  TDOA 
(i.e.,  incoherent  averaging). 

The  performauice  of  the  spectral  domain  based  correlation 
algorithm  developed  in  this  thesis  must  be  further  quantified. 
Sets  of  ROC  curves  should  be  obtained  for  different  SNR'  s 
(i.e.,  SNR  variations  between  channels). 

The  simulated  FFT  error  graph  should  be  validated  by 
measuring  the  error  performance  of  a  dedicated  properly 
dimensioned  FFT. 
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APPENDIX  A.  CORBBLATIOM  MEAN  AND  VARIANCE 


Two  zero  mean,  independent  time  series  are  the  inputs  to 
a  correlator.  This  appendix  will  derive  the  expected  value  amd 
variance  of  the  crosscorrelation  function  output.  In  this 
thesis,  the  noise  is  zero  mean  (i.e.,  shifted)  Rayleigh  noise. 

The  crosscorrelation  function  r^(l)  is  defined  in  Equation 
8  and  is  repeated  here  for  convenience 

H  x(i)y(i  +  J)  (€4) 

i-O 


A.  EXPECTED  VALUE 

Xi  and  are  assumed  to  be  independent,  zero  mean 
sequences 


Elx^]  =  Eiy^}  =  0 . 

The  expected  value  of  r^(l)  is 

w-i-lil  w-i-UI 

E[r^{l)}  =  ^  ^ 

w-i-lil 

=  j;  E{Xi}E[yi.,] 

1-0 

=  0  . 


(65) 


(66) 
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B.  VBRIANCB 


The  random  variable  z  has  a  variance  defined  to  be 

al  =  El(z  -  E{z})^} 

because  E{z}  =0  for  a  zero  mean  sequence ,  (®7) 

al  =  E{z^} . 

The  variance  of  the  crosscorrelation  function  r^(l)  is  then 

w-i-UI  w-i-UI  w-i-UI 

^  E{riy}  =^{(5^  = -PI  E  T. 

i=0  i=0  1=0 

w-i-UI  w-i-ui 

=  E{  T  Y.  . 

j=o 

x^Xj  and  are  two  groups  that  are  independent  of  each 

other.  The  expectation  operation  can  be  applied  to  each  group 
individually . 

w-i-UI  w-i-UI 

"  E  E  EiXiXj)  .  (69) 

1=0 

The  two  indices  i  and  j  define  a  matrix  of  values.  Two 
contributions  have  to  be  considered. 

Contribution  1.  When  ia>j  the  summation  occurs  along  the 
diagonal  of  the  matrix.  Only  one  sviromation  is  used  and 
essentially  the  terms  are  being  squared.  These  squared  terms 
are  not  independent.  To  sixnplify  the  computation  define  m 
equal  to  i  and  j 
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(70) 


N-l-UI 

q  =  E{xl)E{yl,} 


Because  x„  and  y.  are  zero  mean  sequences, 


EKxl)  =  ol. 


and  ai  =  ai  = 


Substituting  <T^  into  Equation  70 


w-i-Ui  w-i-Ul 

Cl  =  E  OxOy  =  5^  1 


=  oM(w-i-)il)  +1] 

=  a*{N-\l\)  . 


Contribution  2.  When  i^j  all  the  terms  except  those  on  the 
diagonal  are  summed.  These  terms  are  independent.  For  i^j 

w-i-Ui  w-i-|i| 

Cj  =  E  E  EixJ  Eix^)  E{yi^^}  E{yJ^^} 

j^o  (72) 

=  0,  using  Equation  65 


Therefore,  Equation  69  becomes  al  =  q  +  q 


=  <J*(N  -  \1\)  . 
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C .  CROSSCOBBBIATIOM  STATISTICS 


Each  data  point  output  by  the  correlation  function 
represents  a  summation  of  terms.  Using  the  central  limit 
theorem,  the  output  of  the  correlation  function  i '  assumed  to 
have  a  Gaussian  distribution. 

In  conclusion  r^(l)  ~  NiO,  o‘(W-|i|) 
where  N  =  Normal  (Gaussian)  distribution  . 
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APPENDIX  B.  MATLAB  CODE 


A.  USERS  GUIDE 

The  frequency  domain  based  correlation  software  given  in 
this  appendix  is  written  in  Pro-MatladD  (version  3.5h).  All 
simulations  were  run  on  a  Naval  Postgraduate  School  Sun  work 
station  using  a  UNIX  operating  system. 


B.  CORRELATION  MATLAB  CODE 

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% 
%  Program  N2une  :  FreqCorr7.m  12  May  91 

% 

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% 

clear 

clg 

clc 

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% 

% 

%  PULSE  CORRELATION 

%  Version  7 . 0 

% 

%  1 .  1  Channel  vs  Reference  Channel 
%  or 

%  1  Cheuuiel  vs  Second  Channel. 

% 

%  2 .  Rician  Signal  and  Rayleigh  Noise  pdf' s  . 

%  (Low  SNR  signals  have  power  approximately  2.) 

% 

%  3.  0  dc  Correlation.  Subtract  dc  from  both  signals. 

% 

%  4.  Normalize  correlation  output  in  Fourier  domain. 
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% 
%  This  MatLah  program  will  either  : 

%  1.  Correlate  a  received  pulse  sequence  against  a  reference 
%  pulse  sequence.  The  reference  sequence  parameters  are 
%  specified  by  the  user  and  are  used  to  construct  the 
%  sequence . 

%  or 

% 
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%  2 .  Correlate  two  received  pulse  sequences  against  each 
%  other. 

% 

% 

%  Correlation  of  the  pulse  sequences  is  accomplished  by 
%  multipling  in  the  frequency  domain  one  spectrum  against 
%  the  conjugate  of  the  other  spectrum.  An  inverse  FFT  is 
%  performed  on  the  result  to  yield  the  time  difference 
%  (TDOA)  between  the  two  sequences. 
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% 


clc 

echo  on 

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% 
%  SELECT  CORRELATION  OPTION 

% 

%  Select  either  1  or  2  : 

% 

%  1 .  Time  delay  in  one  receiver  channel  measured  against  a 

%  reference  signal. 

% 

%  2.  Differential  time  delay  between  two  channels. 

% 

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% 
echo  off 

option  *  input (' Select  correlation  option  1  or  2  ;  '); 

clc 


echo  on 

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% 
%  Option  1  :  Time  Delay  In  One  Chaimel 

% 

%  Reference  Pulse  Train  Pareuneter  Selection. 

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% 
echo  off 
if  option—l 

wi>-input  (  'Pulse  Width  (seconds)  :  '); 

pr«lnput (  'Pulse  Period  (T  seconds)  :  '); 

nu^input (  'Number  of  pulses  in  pulse  train  :  '); 

NuinberPolntsRef«pr*nu; 
end 
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clc 

echo  on 

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% 

% 

% 

%  Delayed  Pulse  Train  Parameter  Selection 

% 

% 

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% 
echo  off 

dl ys input (' Channel  1  delay  tau  (sec.)  of  pulse  train  :  '); 

wi2=input ('Pulse  Width  (seconds)  :  '); 

pr2=input ('Pulse  Period  (T  seconds)  :  '); 

nu2 = input  ('Nvunber  of  pulses  in  Delayed  Pulse  Train  :  ')/ 

NumberPointsDel* (pr2*nu2) +dly;  %calc.  nmbr  pts  sig  +  delay 

clc 

if  option==2 

dl yCh2s input (' Channel  2  delay  tau  (sec.)  pulse  train  :  '); 

wiCh2=input (' Pulse  Width  (seconds)  :  '); 

prCh2=input ('Pulse  Period  (T  seconds)  :  '); 

nuCh2»:input  (' Number  of  pulses  in  Delayed  Pulse  Train  :  '); 

NumberPointsDel2*  (prCh2*nuCh2) +dlyCh2;  %ninbr  pts  sig  +  delay 
end 


%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% 
%  Zero  pad  both  sequences  to  N3sN2+Nl-l.  Must  then  make  N3 
%  meet  the  simple  formula  N3  «  2''m  to  allow  speedy 
%  confutation  of  the  FFT. 

%  WARNING  :  IBM  80286  MATLAB  WILL  NOT  ALLOW  VECTORS  GREATER 
%  THAN  4000.  So  an  input  pulse  train  of  128  will  exceed  the 
%  system.  Solution  is  to  use  the  SUM  workstation  or  386/486. 
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% 
if  option“l 

ZeroPadPoints^NufflberPointsRef+NumberPointsDel'-l;  %M3»N2-i-Nl-l 
for  m«2:l:19;  %2''(19)  =  524288  point  FFT  max 

if  2''m  >■  ZeroPadPoints, 

ZeroPadPoints*2^m; 

break;  %when  find  a  2'^m  power,  exit  loop 

end 

end 

end 
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if  option==2 

ZeroPadPoints=NximberPointsDel2+NuinberPointsDel-l ;  %N3=N2+N1-1 
for  m=2:l:19;  %2'' (19)  =  524288  point  FFT  max 

if  2''m  >=  ZeroPadPoints, 

ZeroPadPoints=2''m; 

break;  %when  find  a  2''m  power,  exit  loop 

end 

end 

end 


%%%%%%%%%%%%%%%  Create  Reference  Pulse  Train  %%%%%%%%%%%%%% 
if  option==l 

ReferencePulse=zeros (1 : ZeroPadPoints) ;  %zero  pad  past  pr*nu 
for  PulseCounter^sl  :pr  :N\afflberPointsRef  %steps  of  Ref.  period 
for  CutUpPulse=0 : pr ; 
if  CutUpPulse<=wi 

ReferencePulse (PulseCounter+CutUpPulse) =1.0; 

end 

end 

end 

Refdc=mean  (ReferencePulse  (1  :NuinberPointsRef) )  ; 

ReferencePulse (1 :NumberPointsRef ) =ReferencePulse (1 rNumberPoi 

ntsRef ) -Refdc; 

end 


%%%%%%%  Create  Chamnel  1  Delayed  Pulse  Train  %%%%%%%%%%%%%% 
DelayedPulseTrain=zeros (1 : ZeroPadPoints) ; %0  pad  past  pr2*nu2 
for  PulseCounter- (dly+1) :pr2: (pr2*nu2+ (dly+l) ) ; %pr2*nu2=Nmbr 
for  CutUpPulse=0 :pr2;  %...pts  in  Del.  Pulse  Train 

if  CutUpPulse<=wi2 

DelayedPulseTrain  (PulseCoiinter-fCutUpPulse)  =1.0; 
end 

end 

end 


if  option==2 

%%%%%%%  Create  Chemnel  2  Delayed  Pulse  Train  %%%%%%%%%%%%%% 
DelayedPulseTrain2=zeros (1 : ZeroPadPoints) ;%0  padpast  pr2*nu2 
for  PulseCounter=  (dlyCh2-t-l)  :prCh2:  (prCh2*nuCh2+ (dlyCh2+l) )  ; 
for  CutUpPulse=0 ; prCh2 ; 
if  CutUpPulse<“wiCh2 

DelayedPulseTrain2 (PulseCounter+CutUpPulse) =1.0; 
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end 


end 

end 

end 


clc 

echo  on 

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% 
%  Signal  to  Noise  Ratio  Calculation 

%  Define  Rayleigh  noise  power  to  equal  2 . 

%  Rician  pdf  created  for  sinusoid  +  Gaussian  noise. 

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% 
echo  off 

SNR=input('  Desired  SNR  ;  '); 

NoisePwr=2;  %rayleigh  noise  pwr  =2 

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% 
%  Sum  Channel  1  &  2  signals  plus  their  initial  delay. 

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% 
SigEnergy»0;  %initialize  Channel  1  energy  to  0 

for  index=l;NumberPointsDel;%sum  over  Ch.  1  delayed  signal 
SigEnergy»SigEnergy+  ( (DelayedPulseTrain  (index) )  .  ''2) ; 

end 

if  option=«2 

SigEnergy2«0 ;  %initialize  Channel  2  energy  to  0 

for  indexsl :yuiiiberPointsDel2; %sum  over  Ch.  2  delayed  signal 
SigEnergy2«SigEnergy2+  ( (DelayedPulseTrain2  (index) )  .''2)  ; 

end 

end 


%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% 

%  Calc.  Ch.  1  signal  amplitude  to  meet  required  input  SNR. 

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% 

SNRIO-SNR/IO; 

SNRLlnear">10^SNR10;  %SNR  in  linear  units 

SigAmplitude^sqrt (NoisePwr*SNRLinear*NumberPointsDel/SigEner 

gy)  ; 


if  option«2 

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% 
%  Calc.  Ch.2  signal  amplitude  to  meet  required  input  SNR. 
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%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% 

SNR10=SNR/10; 

SNRLinear=10^SNR10 ;  %SNR  in  linear  units 

SigAinplitucie2=sqrt  (NoisePwr*SNRLinear*NumberPointsDel2/SigEn 
ergy2)  ; 
end 


%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% 
%  Create  Rayleigh  Noise  matrix  for  Channel  1. 
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% 
rand (' normal' ) /  %Gaussian  mean=0,  variance=l 

for  index=l iNumberPointsDel, 

GaussNoisel (index) =rand;  %N(0,1) 

GaussNoise2 (index) =rand;  %N(0^1) 

end 

RayleighNoise=aqrt  (  (GaussNoisel)  .  ''2+  (GaussNoise2)  .  ''2)  ; 


%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% 
%  Create  Rayleigh  Noise  matrix  for  Channel  2 . 
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% 
if  option=*2 

rand (' normal' ) ;  %Gaussian  mean=0,  variance=l 

for  index»>l :  NmnberPoint sDel2 , 

GaussNoise3 (index) -rand;  %N (0, 1) 

GaussNoise4  ^ndex)  -rand;  %N  (0, 1) 

end 

RayleighNoise2asqrt  ( (GaussNoise3)  .  ^2+  (GaussNoise4)  .  ‘'2) ; 
end 


%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% 
%  Make  Rician  eunplitude  of  Channel  1  delayed  pulse  train 
%  and  remove  the  dc  component  of  entire  signal. 
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% 
SE*0  ; 

for  index>l  zNvimberPointsOel,  %1  — >  delay  +  delayed  signal 
if  DelayedPulseTrain  (index)  »1  %if  one,  make  Rician 


DelayedPulseTrain (index) =DelayedPul3eTrain (index) . *rician (Si 
gAmplitude) ; 

SE= (DelayedPulseTrain (index) ) . ^2+SE; 
end  %end  Rician  modification  loop 
end 
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SP=SE/NuinberPointsDel; 


%%%%%%%  Add  noise  to  Channel  1  Delayed  Pulse  Train  %%%%%%% 
for  index=l iNumberPointsDel,  %1  -->  delay  +  delayed  signal 
if  DelayedPulseTrain (index) =*0  %no  Rayleigh  noise  to  signal 

DelayedPulseTrain (index) =DelayedPulseTrain (index) +RayleighNo 
ise (index) ; 

end 

end 

%%%%%  Remove  dc  component  of  Channel  1  Pulse  Train  %%%% 
Chldc=^ean  (DelayedPulseTrain  (1  :N\imberPointsDel) )  ; 
DelayedPulseTrain (1 :NumberPointsDel) -DelayedPulseTrain (1 :Num 
berPointsDel) -Chide; 


if  option*=2 

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% 

%  Make  Rician  aunplitude  of  Channel  2  delayed  pulse  train 
%  and  remove  dc  component  of  entire  signal. 
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% 
SE«0  ; 

for  index=l :NumberPointsDel2,  %1  — >  delay  +  delayed  signal 
if  DelayedPulseTrain2  (index)  »1  %if  one,  make  Rician 

DelayedPulseTrain2  (index)  sDel’^<yedPulseTrain2  (index)  .  *rician  ( 
SigAmplitude2) ; 

SE®  (DelayedPulseTrain2  (index) )  .  ''2+SE; 
end  %end  Rician  modification  loop 
end 

SP-SE/NumberPointsDel2 ; 


%%%%%%  Add  noise  to  Channel  2  Delayed  Pulse  Train  %%%%%%%% 
for  index«l  :NuiiiberPointsDel2,  %1  — >  delay  +  delayed  signal 
if  DelayedPul8eTrain2  (index)  »0  %no  Rayleigh  noise  to  sig 

DelayedPul8eTraln2 (index) sDelayedPulseTrain2 (index) tRayleigh 
Noise2 (index) ; 

end 

end 

%%%%%  Remove  dc  component  of  Channel  2  Pulse  Train  %%%% 
Ch2dc»mean (DelayedPulseTrain2 (1 :NumberPointsDel2) ) ; 
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DelayedPulseTraln2 (1 :NumberPointsDel2) =DelayedPulseTrain2  (1 : 

NumberPointsDel2)  -Ch2cic; 

end 

subplot  (221)  %two  rows,  two  coliomns 

time= (0 : 1 : ZeroPadPoints-1) ;  %Start  time  axis  at  0 

if  option==l 

%%%%%%%%%%%%%  Plot  Reference  Pulse  Train  %%%%%%%%%%%%%%%%%% 

plot (time, ReferencePulse) ; 

title (' Reference  Pulse  Sequence'); 

xlabel ('time' ) ; 

grid 

end 


%%%%  Plot  Channel  1  pulse  train  +  Rayleigh/Rician  noise%%%% 

topDel=l . l*max (DelayedPulseTrain) ; 

bottomDel«l . l*min (DelayedPulseTrain) ; 

axis([0  NumberPointsDel  bottomDel  topDel]) 

plot (time, DelayedPulseTrain) ; 

title  ('Channel  1  Pulse  Sequence  -f  Rician/Rayleigh  Noise'); 

xlabel  ('time' )  ; 

grid 


if  option==2 

%%%%  Plot  Channel  2  pulse  train  +  Rayleigh/Rician  noise  %%% 
topDel2Bl . l*max (DelayedPulseTrain2) ;  %10  percent  headroom 
bottomDel2Bl . l*min (DelayedPulseTrain2) ; 
axis([0  NumberPointsDel2  bottomDel2  topDel2]) 
plot (time,DelayedPulseTrain2) ; 

title (' Channel  2  Pulse  Sequence  +  Rician/Rayleigh  Noise'); 

xlabel ('time' ) ; 

grid 

end 


%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% 
%  Correlation 

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% 

if  option«l 

fR*fft (ReferencePulse) ; 

fD»fft (DelayedPulseTrain) ; 

end 
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if  option==2 

fR=fft (DelayedPulseTrain) ; 
fD=fft (DelayedPulseTrain2) ; 
end 

z=fD . *con j (fR) ; 

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% 
%  Calculate  normalizing  coefficient 
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% 
Norml=0; 

Norm2=0; 

for  k=l : 1 : length (fD) 

Norrnl®  (aba  (fD  (k) )  .  ''2)  +Norml; 

Norm2=  (a±>s  (fR  (k/  )  .  ''2)  +Norm2; 

end 

Norml=sqrt (Morml) ; 

Norm2*sqrt (Norm2) ; 

NormCoef f=Norml *Norm2; 
z*z . /NormCoef f / 

z=z*length(fD) ; 

ifftz=ifft (z) ; 
magz»abs (ifftz) ; 


%%%%%  Flip  vector  for  appearance  %%%%%%%%%%%%%%% 
Topmagza (length (magz) /2+1) ; 

Bottommagzslength (magz) /2; 

magzLen«length (magz) ;  %measure  length  of  magz 
Transmagz* [magz (Topmagz :magzLen) ,magz (1 :Bottommagz) ] ; 

%%%%%%  Make  -time  to  +time  axis  %%%%%%%%%%%%%%%% 
minuatime"- (length (Transmagz) /2) ; 
posittime" (length (Transmagz) /2) -1; 
timel> (minustime: 1 iposittime) ; 

topmagz«l . l*max (magz) ;  %add  10  percent  head  room 
axis ( [minustime  posittime  0  topmagz])  %scale  axis 

subplot (212) 

plot (timel, Transmagz) ; 

title ('FFT  Correlation  :  Channel  1  vs.  Channel  2'); 


79 


xlabel ( ' time  lag'); 
grid 

gtext (' +7dB  SNR' ) 

lj4print  %Spanagel  Room  427  Laserjet 

axis<[l  2  3  4]);  %reset  axis  scaling 
axis; 


%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% 
%  Prograun  Naune  :  CorrPfaLoop.m  19  July  91 

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% 
%  Purpose  :  To  compute  the  Pfa  for  Signal  /  Noise  only  for 
%  then  correlation  function.  This  program  creates  zero  mean 
%  RAYLEIGH  noise. 

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% 

clear 

clg 

subplot (221) 
raund(' normal' ) 

PFA»1.0;  %will  be  divided  by  ten  to  start 

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% 
for  PFACOUNTER-1 ; 1 ; 4 
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% 

PFA«PFA/10;  %Pfa  -  0.1  0.01,  0.001,  0.001 

%%%%%%%%%%%%%%%%%%%%%%%%%%%% 

for  SNR«0:1:18  %Outer  loops  through  all  SNR's. 

%%%%%%%%%%%%%%%%%%%%%%%%%%%% 

Iterations>"30;  %30  data  points  per  SNR  point 

%%%%%%%%%%%%%%%%%%%%%%%%%%% 

for  loop«l : 1 : Iterations  %inner  loop,  creates  groups  of  +-  a 

%%%%%%%%%%%%%%%%%%%%%%%%%%% 

loop; 

AboveThreshold(loop) «0; 

BelowThreshold ( loop) ^O ; 


Figur*  24.  FreqCorrV.m  MATLAB  flowchart 


Threshold (loop) =0; 

MaxRxy (loop) =0; 

RxyThreelag (loop) =0; 

%%%%%%%%%%%%  Set  Parameters  %%%%%%%%%%%%%%%%%%%%%%%% 

N=100;  %100  noise  points,  therefore  200  correlation  points 

a=10;  %+-  10  %  on  either  side  of  0  time  lag 

Pfa=PFA;  %ReestadDlish  because  Pfa  destroyed  each  loop 
PfaDef ined=Pfa; 

%%%%%%%%%%%%%  Make  N  points  Noise  %%%%%%%%%%%%%%%%%% 
snrnoiseloop  %call  signal  plus  noise  generation  at  a  SNR 


%%%%%  Correlation  of  zero  mean  noise  only  for  Pfa  calcul.  %% 
MeanRayl^mean (RayleighNoise) ;  %noise  created  in  snrnoise.m 
MeanRay2>mean (RayleighNoise2) ; 
RayleighNoise^RayleighNoise-MeanRayl; 
RayleighNoise2-RayleighNoise2-MeanRay2; 

Rxyaxcorr (RayleighNoise, RayleighNoise2) ;  %xcorr  signal+noise 
later 

%%%%%%%%%%%%%%%%%%  Compute  t"  %%%%%%%%%%%%%%%%%%%%%%% 
t*0; 

for  i=(N-a) !l: (N+a) 
t»t+Rxy (i) ; 

end 

t=t/ (2*a+l) ; 

%%%%%%%%%%%%%%%  Compute  variemce  t  %%%%%%%%%%%%%%% 

Vart-0 ; 

for  i= (N-a) ; 1 : (N+a) 

Vart»Vart+  (Rxy  (i)  -t)  .  ''2; 

end 

Vart-Vart/ (2*a+l) ; 

%%%%%%%%%%%%%%%%%  Define  Pfa  %%%%%%%%%%%%%%%%%%%% 

Pfa*Pfa*2;  %multiply  by  2 

Pfa»l-Pfa;  %svibtract  one 

Threshold (loop) "inverf (Pfa) ; 

Threshold (loop) -Threshold (loop) *sqrt (2) *sqrt (Vart) ;  %mult  by 
the  sqrt (2) 

%%%%%  Xcorr  of  signal  +  zero  mean  noise  %%%%%%%%%% 

Rxy-xcorr (DelayedPulseTrain,DelayedPulseTrain2) ; 
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%%%%%%%  Compare  +a  to  -a  to  Threshold  %%%%%%%%%%%% 
for  i= (N-a) : 1 : (N+a) 

if  Rxy (i) >=Threshold (loop) 

AboveThreshold (loop) =AboveThreshold (loop) +1; 
else 

BelowThreshold (loop) «BelowThreshold (loop) +1/ 

end 

if  Rxy (i) >MaxRxy (loop) 

MaxRxy  (loop)  =“Rxy  (i)  / 

SaveThreshold (loop) “Threshold (loop) ; 

%SaveVart (loop) ®sqrt (Vart) ; 

end 

end 

RxyThreelag (loop) =Rxy (97) ; 

end  %end  #  Iterations  loop 

%%%%%%%  Count  correct  threshold  crossings  %%%%%%%% 
Correct TDOA=0 ; 
for  count=l ; 1 : Iterations 

i f  Rxy Three 1 ag ( count ) “*MaxRxy ( count ) 

Correct  TDOA«Correct  TDOA+ 1 ; 
end 

end 

PercentCorrectTDOA* (CorrectTDOA*100) /iterations; 

diary  flag4 
PercentCorrectTDOA 
BelowThreshold 
AboveThreshold 
PfaiDefined 
SNR 

S  aveThr e  sho 1 d 
lag  below 
RxyThreelag 
three  pulses. 

MaxRxy 
diary  off 

end  %end  #  test  points  loop 

end  %end  Pfa  loop 


%compare  threshold  to  zeroth 
%Channel  2  lags  Channel  1  by 
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%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% 
%  Program  Neune  :  snrnoiseloop.m  19  JULY  91 

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% 
dly=0;  %Channel  1  delay 

wi2=5;  %Channel  pulse  width 

pr2=10;  %Channel  1  pulse  period 

nu2=10;  %Channel  1  ntimber  of  pulses 

NumberPointsDel= (pr2*nu2) +dly;  %calc.  nmbr  pts  sig  +  delay 

dlyCh2=3;  %Channel  2  delay 

wich2=5;  %Channel  2  pulse  width 

prCh2-10;  %Channel  2  pulse  period 

nuCh2=10;  %Channel  2  number  of  pulses 

NumberPointsDel2“ (prCh2*nuCh2) +dlyCh2;  %nmbr  pts  sig  +  delay 

ZeroPadPoints=NuinberPointaDel2+NumberPointaDel-l;  %N3=N2+N1-1 
for  m=2:l:19;  %2'‘(19)  =  524288  point  FFT  max 

if  2''m  >*  ZeroPadPoints, 

ZeroPadPoints»2^m; 

break;  %when  find  a  2^m  power,  exit  loop 

end 

end 

%%%%  Create  Channel  1  Delayed  Pulse  Train  %%%%%%%%%%% 
DelayedPulaeTrain=zeros (1 : ZeroPadPoints) ; %0  pad  past  pr2*nu2 
for  PulseCounter»(dly+l) :pr2; (pr2*nu2+ (dly+1) ) / %pr2*nu2=Nmbr 
for  CutUpPulse»0:pr2;  %...pts  in  Del.  Pulse  Train 

if  CutUpPulse<»wi2 

DelayedPulseTrain (PulseCounter+CutUpPulse) =1.0; 

end 

end 

end 

%%%%%  Create  Channel  2  Delayed  Pulse  Train  %%%%%%%%%% 
DelayedPulseTrain2=zeros (1: ZeroPadPoints) ;%0  padpast  pr2*nu2 
for  PulseCounter*  (dlyCh2-»-l)  ;prCh2;  (prCh2*nuCh2+ (dlyCh2+l) )  ; 
for  CutUpPulse=0 : prCh2 ; 
if  CutUpPul8e<=wiCh2 

DelayedPulseTrain2 (PulseCounter+CutUpPulse) =1.0; 
end 

end 

end 


%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% 
%  Signal  to  Noise  Ratio  Calculation 

%  Define  Rayleigh  noise  power  to  equal  2 . 


%  Rician  pdf  created  for  sinusoid  +  Gaussian  noise. 

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% 
echo  off 

NoisePwr=2;  %rayleigh  noise  pwr  =2 

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% 
%  Sum  Channel  1  &  2  signals  plus  their  initial  delay. 

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% 
SigEnergy=0;  %initialize  Channel  1  energy  to  0 

for  index=l tNumberPointsDel; %sum  over  Ch.  1  delayed  signal 
SigEnergy=SigEnergy+  ( (DelayedPulseTrain  (index) )  .  ''2)  ; 

end 

SigEnergy2=0;  %initialize  Channel  2  energy  to  0 

for  index=l :NumberPointsDel2; %sum  over  Ch.  2  delayed  signal 
SigEnergy2=SigEnergy2+ ( (DelayedPulseTrain2 (index) ) . ^2) ; 

end 

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% 

%  Calc.  Ch.  1  signal  aunplitude  to  meet  required  input  SNR. 
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% 
SNR10»SNR/10; 

SNRLinear«10'^SNR10;  %SNR  in  linear  units 

SigAmplitude«sqrt (NoisePwr*SNRLinear*NumberPointsDel/SigEner 

gy) ; 

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% 

%  Calc.  Ch.2  signal  amplitude  to  meet  required  input  SNR. 

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% 

SNRlO-SNR/10; 

SNRLinear^lO'^SNRlO;  %SNR  in  linear  units 

SigAmplitude2asqrt (NoisePwr*SNRLinear*NumberPointsDel2/SigEn 
ergy2) / 

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% 

%  Create  Rayleigh  Noise  matrix  for  Channel  1. 
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% 
for  index^l :NumberPointsDel, 

GaussMoisel (index) «rand;  %N(0,1) 

GaussNoise2 (index) «rand;  %N(0,1) 

end 

RayleighNoise^sqrt  ( (GaussNoisel)  .  ^2-f  (GaussNoise2)  .  ''2)  ; 

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% 

%  Create  Rayleigh  Noise  matrix  for  Channel  2 . 
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% 
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%N(0, 1) 
%N(0, 1) 


for  index=l  :NuinberPointsDel2, 

GaussNoise3 (index) =rand; 

GaussNoise4 (index) =rand; 

end 

RayleighNoise2=3qrt  (  (GaussNoiseS)  .  ''2+  (GaussNoise4)  .''2)  ; 

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% 

%  Make  Rician  amplitude  of  Channel  1  delayed  pulse  train 
%  and  remove  the  dc  component  of  entire  signal. 
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% 
SE=0; 

for  index=:l  :NumberPointsDel,  %1  -->  delay  +  delayed  signal 
if  DelayedPulseTrain  (index)  »1  %if  one,  make  Rician 

DelayedPulseTrain (index) =DelayedPulseTrain (index) . *rician (Si 
gAmplitude) ; 

SE=:  (DelayedPulseTrain  (index) )  .  ^2+SE; 
end  %end  Rician  modification  loop 

end 

SP=SE/NumberPoint sDel ; 

%%%%%%%  Add  noise  to  Chamnel  1  Delayed  Pulse  Train  %%%%%%% 
for  index-1 ;NumberPointsDel,  %1  — >  delay  +  delayed  signal 
if  DelayedPulseTrain  (index)  “0  %no  Rayleigh  noise  to  signal 

DelayedPulseTrain (index) -DelayedPulseTrain (index) +RayleighNo 
ise (index) ; 

end 

end 

%%%%%  Remove  dc  component  of  Channel  1  Pulse  Train  %%%% 
Chldc-meeui  (DelayedPulseTrain  (l:NvunberPointsDel) ) ; 
DelayedPulseTrain  (1  :NumberPointsDel)  -DelayedPulseTrain  (1  :Nvim 
berPointsDel) -Chide; 

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% 

%  Make  Rician  euaplitude  of  Channel  2  delayed  pulse  train 
%  and  remove  dc  component  of  entire  signal. 
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% 
SE-0; 

for  index-1 :NumberPointsDel2,  %1  — >  delay  +  delayed  signal 
if  DelayedPulseTrain2  (index) —1  %if  one,  make  Rician 

DelayedPulseTrain2 (index) -DelayedPulseTrain2 (index) . *rician ( 
SigAmplitude2) ; 

SE- (DelayedPulseTrain2 (index) ) . ^2tSE; 


end  %end  Rician  modification  loop 
end 

SP=SE/NumberPointsDel2 ; 

%%%%%%  Add  nc  3e  to  Channel  2  Delayed  Pulse  Train  %%%%%%%% 
for  index=l  :N’JimberPointsDel2,  %1  — >  delay  +  delayed  signal 
if  DelayedPulseTrain2 (index) ==0  %no  Rayleigh  noise  to  sig 

DelayedPulseTrain2 (index) =DelayedPulseTrain2 (index) +Rayleigh 
Noise2 (index) ; 

end 

end 

%%%%%  Remove  dc  component  of  Channel  2  Pulse  Train  %%%% 
Ch2dc=mean (DelayedPulseTrain2 (1 :NumberPointsDel2) ) ; 
DelayedFulseTrain2  (1  :NuinberPointsDel2)  =DelayedPulseTrain2  (1 : 
NumberPointsDel2) -Ch2dc; 

DelayedPulseTrain=DelayedPulseTrain (1 : 100) ; 
DelayedPulseTrain2*DelayedPulseTrain2 (1 : 100)  ; 

RayleighNoise^RayleighNoise (1 : 100) ; 
RayleighMoise2»RayleighNoise2 (1 : 100) / 
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Figur«  25.  CorrPfaLoop.m  MATLAB  flowchart. 


APPENDIX  C.  TDOA  CONSTANT  THRESHOLD  SIMULATION 


X.  INTRODUCTION 

Four  graphs  of  the  output  of  a  time  domain  based 
correlation  detector  using  a  constant  threshold  are  given. 
Each  figure  has  four  simulations  using  a  SNR  of  IdB.  The 
is  varied  from  0.01  to  0.00001  from  Figure  26  to  Figure  29. 
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Figur*  26.  Time  domain  based  correlation  detector  output 
for  a  SNR  of  IdB  and  a  P,.  of  0.01. 


time  domain  crosscorrelntion  with  threshold 


Figure  27.  Time  domain  based  correlation  detector  output 
for  a  SNR  of  IdB  and  a  Pj.  of  0.001. 
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time  domain  crosscorrelatiun  with  threshold  time  domain  crosscorrelation  with  threshold 
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time  domain  crosscorrelation  with  threshold  time  domain  crosscorrelation  with  threshold 


(l)^ 


Figure  29.  Time  domain  based  correlation  detector  output 
for  a  SNR  of  IdB  and  a  Pf.  of  0.00001. 
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t  ldB  SNR,  pia=().lM)0|)l 


APPENDIX  D.  FFT  OUTPUT  NOISE  TO  SIGNAL  RATIO 

A.  INTRODUCTION 

The  noise  to  signal  ratio  for  thirteen  passes  through  a 
FFT  butterfly  are  calculated  below.  The  calculation  assumes 
the  word  length  of  the  arithmetic  logic  unit  (ALU)  of  the 
computer  is  16  bits  long.  The  length  of  the  input  data  points 
are  12  bits  long.  The  data  points  are  right  justified  when 
placed  into  the  ALU.  For  a  16  bit  ALU  and  a  12  bit  data 
input,  scaling  occurs  at  the  seventh  pass  and  every  other  pass 
after  that. 

Pass  1 

The  noise  power  at  the  output  of  the  butterfly  after  the 
first  pass  N*^(l)  is  a  function  of 

nI(1)  =2iVi(0)  +  2o|(0)  +  S2ot(0)  +  ot.(0) 

where  =  signal  power  input  to  butterfly 

Op  =  truncation  noise  power  ^74) 

o|  =  scaling  noise  power 
ol,  =  sine/ cosine  table  noise  power 
Ari(O)  =  noise  power  input  to  FFT  . 

During  pass  one  there  is  no  scaling,  no  truncation,  and  the 
trig  tzUble  is  not  used.  These  error  terms  are  zero.  Equation 
74  reduces  to 
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(75) 


N^il)  =  2N~ 

where  =  ^2'^® 

=  noise  input  to  butterfly 


Pass  2 

The  noise  power  at  the  output  of  the  butterfly  after  the 
second  pass  N,^^(2)  is  a  function  of 

wi(2)  =  2Wi(l>  +  2o|{l)  +  252ot(l)  +  Or(l>  .  (76) 

During  pass  two  scaling  and  truncation  are  not  performed,  and 
the  trig  teible  is  not  used.  These  noise  terms  are  zero. 
Substituting  Equation  75  into  76 


iVi(2)  =  2(2iV^) 
=  . 


(77) 


Pass  3 

The  noise  power  at  the  output  of  the  butterfly  after  the 
third  pass  Nj^^O)  is  a  function  of 

N^(3)  =2Ari(2)  -f  2o|(2)  +225^0^(2)  +0^(2)  .  (78) 

Scaling  noise  is  zero  because  it  is  not  performed  during  pass 
three.  Substituting  Equation  77  into  78 
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wiO)  =  2(2^^,)  +  2^S^al  +  <Jt 


(79) 


=  +  225^0",  +  al  . 


Pass  4 

The  noise  power  at  the  output  of  the  butterfly  after  the 
fourth  pass  Nx^(4)  is  a  function  of 

ivi(4)  =2W|(3)  +  2o|(3)  +  2^S^al{3)  +o|(3)  (80) 

Scaling  noise  is  zero  because  it  is  not  performed  during  pass 
four.  Substituting  Equation  79  into  80 

//i(4)  =  2^Ng  +  2*S^al  +  3al  .  (81) 


Pass  5 

The  noise  power  at  the  output  of  the  butterfly  after  the 
fifth  pass  N,^*(5)  is  a  function  of 

N^{5)  =2Ni(4)  +  2o|(4)  *  +0^(4)  .  (82) 

Scaling  noise  is  zero  because  it  is  not  performed  during  pass 
five.  Substituting  Equation  81  into  82 

iy|(5)  =  2^Ng  +  (2’  +  2*)S^al  +  70r  ■  <®3) 


Pass  6 

The  noise  power  at  the  output  of  the  butterfly  after  the 
sixth  pass  Nj^^(6)  is  a  function  of 
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wi(6)  =2Ni(5)  +20^(5)  +2-520^(5)  +  oi(5)  .  (8^) 

Scaling  noise  is  zero  because  it  is  not  performed  during  pass 


six.  Substituting  Eq[uation  83  :-nto  84 

Wj(6)  =  2^Nq  >  2'’S^ai  +  150r  •  <85) 


Pass  7 

The  noise  power  at  the  output  of  the  butterfly  after  the 
seventh  pass  N;^*(7)  is  a  function  of 

A/i(7)  =2Ni(6)  +  2o|(6)  ^■2«S^ot(6)  +0^(6)  .  (86) 

Scaling  is  performed  for  the  first  time  during  pass  seven. 
Sxibstituting  Equation  85  into  86 

nI(7)  *  2’W^  +  (2®  +  2^)S^ol  +  3lOr  +  2o|  .  (87) 


Pass  8 

The  noise  power  at  the  output  of  the  butterfly  after  the 
eighth  pass  Nj^^(8)  is  a  function  of 

NliB)  =  2NI{1)  +  2o|(7)  *2'’S^oiO)  +0^(7)  .  (88) 

Scaling  is  not  performed  during  this  pass.  Only  the  scaling 
noise  from  the  previous  pass  is  added.  Substituting  Equation 
87  into  88 

Wi(8)  =  2^Ng  +  (2®  *  2^)S^oi  +  63o|  +  40^  .  (89) 
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Pass  9 


The  noise  power  at  the  output  of  the  butterfly  after  the 
ninth  pass  (9)  is  a  function  of 

Nil9)  =  2Ni(8)  +20^(8)  --2®S-ot(8)  +oi(8)  (90) 

Substituting  Equation  89  into  90 

Wi(9)  =  2®W^  +  (2^°  +  2®  +  2*)S^ot  +  1270r  +  lOo,  .  OD 


Pass  10 

The  noise  power  at  the  output  of  the  butterfly  after  the 
tenth  pass  Nj^^<10)  is  a  function  of 

wi(lO)  =  2nI(9)  +  20s(9)  +2*520^(9)  +0^(9)  (92) 

Scaling  is  not  performed  this  pass.  Only  the  scaling  noise 
from  the  previous  pass  is  added.  Substituting  Equation  91  into 
92 


Ari(lO)  =  2^°iV^  +  2^252^2  +  2550t  +  20a\  .  (93) 

Equation  93  describes  the  noise  power  at  the  output  of  a 
butterfly  after  the  tenth  data  pass.  The  noise  to  signal  ratio 
after  the  tenth  pass  through  the  butterfly  is  expressed  as  the 
ratio 


( 
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(94) 


Pass  11 

The  noise  power  at  the  output  of  the  butterfly  after  the 
ll“*  pass  Njk*(ll)  is  a  fxinction  of 

W^dl)  •  2Nji(10)  +20^(10)  ♦  2^®S*ai(10)  +0^(10)  (96) 

Substituting  Equation  93  into  96 

Ni(ll)  =  2“nJ  +  (2^^  *  2^'»)S*oi  ♦  511<4  *  42oi  .  (97) 
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Pass  12 


The  noise  power  at  the  output  of  the  butterfly  after  the 
12'"*'  pass  N;^^(12)  is  a  function  of 

Ni(12)  =  2Ni(ll)  +  20s(ll)  +  2^"s"<y»(ll)  *0,(11)  •  08) 

Scaling  is  not  performed  this  pass.  Only  scaling  noise  from 
the  previous  pass  is  added.  Siabstituting  Equation  97  into  98 

Ni(12)  =  2^^N*  *  (2^«  *  2“)S*oi  +  10230^  *  840^  .  (99) 

Pass  13 

The  noise  power  at  the  output  of  the  butterfly  after  the 
IS^**  pass  Nj^^(13)  is  a  function  of 

(100) 

Ni(13)  =  2Wi(12)  +  2<t|(12)  *  2“S*ai(12)  *0^(12)  . 

Substituting  Equation  99  into  100 

(101) 

Ni(13)  =  +  (2«  +  2“  *  2^*)S2ai  +  20470^  +  170c^  . 

Equation  101  describes  the  noise  power  at  the  output  of  a 
butterfly  after  the  IS***  data  pass.  The  noise  to  signal  ratio 
after  the  13^''  pass  is  expressed  as  the  ratio 

2”nJ  +  (2“  +  2“  +  2“)S*ai  +  2047<4  +  170o^  (102) 

•F*  ’  "  - 2^ -  ■ 

Dividing  through  by  the  denominator,  each  independent  and 
uncorrelated  noise  term  can  be  identified 


nI  ^  2  1  I  <^s 

-p  (13)  =-p  "  "  -4-32  *  ■W7I~p 


where  —1(13)  =  output  Noise  to  Signax  ratio 
N^ 

_I  =  input  Noise  to  Signal  ratio 
S* 

5.5aw  =  sine/ cosine  table  noise 
_2 

=  truncation  Noise  to  Signal  ratio 


(103) 


_ f  =  scaling  Noise  to  Signal  ratio  . 

TFT?  S" 


101 


LIST  OF  REFERENCES 


1.  Sonnenberg  G.J.,  Radar  and  Electronic  Navigation.  D.  Van 
Nostrand  Company,  INC.,  1951. 

2.  Adamo  R.C.,  Kalman  Filtering  in  the  Spectral  Domain. 
Master's  Thesis,  Naval  Postgraduate  School,  Monterey, 
California,  March  1991. 

3.  Skolnik  M.I.,  Radar  Handbook.  McGraw-Hill,  1990. 

4.  Skolnik  M. I.,  Introduction  to  Radar  Systems.  McGraw-Hill, 
1980. 

5.  Oppenheim  A. V. ,  Digital  Signal  Processing.  Prentice-Hall, 
1975. 

6.  Bendat  J.S.  and  Piersol  A.G.,  Random  Data  Analysis  and 
Measurement  Procedures .  Wiley-Interscience,  1986. 

7.  harris  f.j..  Working  papers,  2  May  1991. 

8.  harris  f.j.,  Working  papers,  9  August  1991. 


INITIAL  DISTRIBUTION  LIST 


1.  Defense  Technical  Information  Center  2 

Cameron  Station 
Alexandria,  VA  22304-6145 

2.  Library,  code  52 
Naval  Postgraduate  School 
Monterey,  CA  93943-5002 

3.  Chairman,  Code  EC 
Department  of  Electrical  and 
Naval  Postgraduate  School 
Monterey,  CA  93943-5000 

4.  Professor  Ralph  Hippenstiel, 

Department  of  Electrical  and 
Naval  Postgraduate  School 
Monterey,  CA  93943-5000 

5.  Professor  Harold  Titus,  Code 
Department  of  Electrical  and 
Naval  Postgraduate  School 
Monterey,  CA  93943-5000 

6.  Professor  Roberto  Cristi,  Code  EC/Cx  l 

Department  of  Electrical  and  Computer  Engineering 
Naval  Postgraduate  School 

Monterey,  CA  93943-5000 

7.  Professor  Charles  Therrien,  Code  EC/Th  1 

Department  of  Electrical  and  Computer  Engineering 
Naval  Postgraduate  School 

Monterey,  CA  93943-5000 

8.  Mr.  Frank  Mika  2 

3246  Test  Wing 

Electronic  Warfare  Test  Division,  TZW 
Eglin  AFB,  FL  32542-5000 

9.  Mr.  George  Barrow  1 

3246  Test  Wing 

Range  Systems  Directorate,  TFRR 
Eglin  AFB,  FL  32542-5000 


2 


1 

Computer  Engineering 


Code  EC/Hi  3 

Computer  Engineering 


EC/Ti  1 

Computer  Engineering 


103 


1 


10,  Naval  Ocean  Systems  Center 

ATTN:  Dr.  C.E.  Persons,  Code  732 
San  Diego,  CA  92152 


104 


